Ethics statement

All methods were carried out in accordance with relevant guidelines and regulations. All experimental protocols were approved by the Animal Welfare Committee of Agro-biotechnology of China Agricultural University. All animals used in this study were cared for and experimented on according to the requirements of the Animal Welfare Committee of Agro-biotechnology of China Agricultural University with the approval SKLAB-2014-06-07.

Experimental population and sample preparation

We aim to assess ddGBS performance in a population with a low level of LD. A nine-generation advanced intercross population was established from two divergent chicken lines, High Quality chicken Line A (HQLA), a broiler line bred by Guangdong Wiz Agricultural Science and Technology, Co. (Guangzhou, China), and Huiyang Beard chicken (HB), a native Chinese meat-type breed. The F0—F₂ cross has been described in detail by Sheng et al. [31]. After F₂ generation, the population was bred by random mating. In total, a set of animal material, consisting of 31 F₀ individuals, 191 F₈ animals, and 602 F₉ progeny, was selected. DNA was extracted from EDTA-anticoagulated blood using the Qiagen DNeasy Blood and Tissue Kit according to the manufacturer’s instructions (Qiagen, Hilden, Germany).

Pre-sequencing processing and evaluation

In this study, both in silico simulation and empirical evidence were considered in choosing the proper enzyme for the digestion of chicken genome. We employed seven single-enzyme digestions (EcoR I, HinP1 I, ApeK I, Pst I, Mse I, Msp I, and Bgl II, including four-/five-/six-cutter enzymes and restriction enzymes resistant to dam, dcm, and CpG methylation, methylation-sensitive: EcoR I, HinP1 I, and ApeK I; methylation-insensitive: Pst I, Mse I, Msp I, and Bgl II) and eight double-enzyme digestions (Pst I-Mse I, Pst I-ApeK I, EcoR I- Mse I, Bgl II- ApeK I, Pst I- Msp I, HinP1 I- Mse I, HinP1 I- ApeK I, and EcoR I- Msp I) in our experiment. In silico analysis were conducted with an in-house Perl script. The size distribution of enzyme digestion fragments was reported using R software. Enzyme digestion experiments for all enzyme digestion combinations were performed according to the enzyme manufacturer’s protocol (New England Biolabs, Ipswich, MA, U.S.A.), and the digesting time for each combination was either 2 h or 12 h.

We employed double-digest genotyping by sequencing (ddGBS) on three samples from the F₀ generation. All DNA samples were diluted to 50 ng/μL, and 200 ng DNA was used for each digestion of the eight double-enzyme combinations according to the enzyme manufacturer’s instructions. We designed 24 barcode adapters (eight enzyme combinations × three samples, see S1 Table). Mixing proportions of the barcode adaptors (BAs) and common adapters (CAs) were determined according to the fragment counts resulted from the in silico analysis of each restriction enzyme combination (S1 Protocol). The barcode adaptors (BA) were linked to the reverse complementary sequences of the Enzyme I overhang, and the common adaptors (CA) were linked to the reverse complementary sequences of the Enzyme II overhang. Library size-selection was implemented by Agencourt^® AMPure^® XP Reagent (Beckman Coulter, Pasadena, CA, U.S.A.): 0.8× and 1.3× sample volume of Agencourt^® AMPure^® XP Reagent can remove most of the short fragments (< 300 bp) and long fragments (> 650 bp), respectively. Detailed library preparation procedures are provided in S1 Protocol.

We evaluated each double-enzyme digestion strategy based on the enzyme digestion fragment size, the fragment consistency index (FCI), the coefficient of variation of sequencing depth (per fragment) across three samples (CV_depth), the number of SNPs, and the distribution uniformity of SNPs across the chromosomes. We also subsampled reads of each individual in different proportions (10%, 20%, 50%, 80%, and 100%), and evaluated the “optimal read-count point (ORP)” and the “saturated read-count point (SRP)” parameters for cost optimization. A detailed definition of the above technical terms was described in the “Terminology” section.

EcoR I- Mse I library preparation

All DNA concentrations were normalized to 50 ng/μL. Samples were digested for 12 h at 37°C with EcoR I- Mse I (New England Biolabs, Ipswich, MA, U.S.A.) in 20 μL volume containing 4 μL DNA (200 ng), 1× CutSmart^® Buffer, 5U EcoR I, and 5U Mse I. The enzymes were then inactivated by heating at 65°C for 20 min, and the samples were cooled to 4°C. The barcode adaptor (EcoR-BA) binds to the EcoR I overhang, and the common adaptor (Mse-CA) matches the Mse I overhang. The 96 indexes at the 3’ end of the barcode adaptors were designed by the GBS Barcode Generator (http://www.deenabio.com/) and modified to allow for Illumina NextSeq500 sequencer (San Diego, CA, U.S.A.) (no barcodes begin with GG; S2 Table). Barcodes were modulated in length between six and nine bases to prevent a decrease of sequencing quality near the restriction sites. 5 μL anneal adapter mix (the ratio of the EcoR-BA and Mse-CA is 0.8:15 based on the predicted fragment counts obtained from EcoR I and Mse I, S3 Table) was ligated to 20 μL digestion products by T4 DNA ligase (Invitrogen, Carlsbad, CA, U.S.A.). The reaction was incubated at 22°C for 1 h, and inactivated at 65°C for 20 min. Considering the maximum reads per flow cell of the NextSeq500 sequencer and the ORP of EcoR I- Mse I, 96 ligation products were pooled together (one library). Agencourt^® AMPure^® beads (Beckman Coulter, Pasadena, CA, U.S.A.) were used for DNA fragment purification and size-selection. The PCR amplification reaction system contained 10 ng purified products, 50 μL Platinum^® PCR SuperMix High Fidelity (Thermo, MA, U.S.A.), and 25 pmol primers (S1 Table). The amplification cycling protocol was as follows: 95°C for 5 min; three steps of 95°C for 30 s, 62°C for 30 s, and 68°C for 30 s for 17 cycles; followed by a final extension at 72°C for 5 min. PCR products were also purified by Agencourt^® AMPure^® beads. The fragment sizes obtained by this method were approximately 300 bp-650 bp, and the fragment size of the highest proportion was 350 bp. The final library quality (concentration and fragment size distribution) was determined by Qubit2.0 Fluorometer (Thermo, MA, U.S.A.) and Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, U.S.A.), respectively.

Sequencing and data processing

All sequencing experiments were performed on the Illumina Nextseq500 Sequencer at the State Key Laboratory for Agro-biotechnology, China Agricultural University. BCL files as primary sequencing output were converted into FASTQ files using bcl2fastq2 conversion software (version 2.16.0). During the conversion step, we also masked and trimmed the sequencing adapter [32]. After the trimming step, the Illumina 91-bp single-end reads were subjected to a filtering process: at first, the reads that were polluted by the adapter sequence were deleted, and then the reads which contained more than 50% low quality bases or more than 5% N bases were removed. The quality control check report of filtered reads was generated by FastQC software (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). We used TASSEL GBS analysis pipeline (version 4.0) [11,33], in which reads were aligned to the chicken reference genome Gallus_gallus-4.0 (released 2011) using Bowtie2 [34]. All SNP filter options in TASSEL were "-c 3", the minimum number of times a tag must be present to be output; "-mnTCov 0.01", the minimum SNP call rate for a taxon to be included in the output; "-mnSCov 0.6", the minimum sample call rate for a SNP to be included in the output; and "-mnMAF 0.05", the minimum minor allele frequency. The raw SNP sites were filtered by VCFtools [35] according to the following parameters: 1) minor allele frequency (MAF) > 5%; 2) genotypes with a quality above 98 (GQ ≥ 98) and depth ≥ 5; 3) and only biallelic markers were retained. Ungenotyped markers were imputed using Beagle4.0 software [36] with the pedigree file of F₈-F₉ family relationships. To annotate mutations from the GBS output, we used the SNPEff program [37], with the chicken reference genome sequence and GTF annotation files downloaded from Ensembl (http://www.ensembl.org/info/data/ftp/index.html). The Circos software package (http://circos.ca/) [38] was utilized to visualize the distribution of fragments, GC islands, repeat regions, and SNPs in the chicken genome. The genome-wide LD pattern assessment was implemented using a squared allelic correlation coefficient (r²) against the distance between the SNPs. To visualize the LD pattern, the r² values were plotted against the pair-wise SNP distances.

Terminology

A “good barcode read” is a sequence read with a perfect match to one of the barcodes provided in a barcode file. A “tag” refers to a unique sequence (excluding the barcode) from one or more “good barcode reads”. A “fragment” is defined as a set of tags that align to the exact same genomic position and strand. The number of tags and fragments is counted by the output file of the TASSEL software [33]. The fragment consistency index (FCI) is defined as the average fragment count from three samples divided by the total fragment counts obtained from pools of three samples. The sequencing depth is calculated as the total good barcode read counts divided by the fragment counts. The CV_depth is calculated as the mean of sequencing depth (per fragment) across three samples divided by the standard deviation (SD). The SNP density is defined as SNP number divided by chromosome length. The CV_{SNP density} is calculated as the mean of SNP density (per chromosome) divided by the standard deviation (SD).

The sequencing cost per fragment unit is calculated by total sequencing cost against fragment counts. The optimal read-count point (ORP) is defined as the minimum sequencing cost per fragment unit. We also defined the saturated read-count point (SRP) as the minimum good barcode reads when reaching the maximum fragment counts.

Screen the appropriate enzyme combinations for the chicken genome

We performed a series of assessments for enzyme selection. The first parameter tested was fragment size. According to the predicted results from the in silico digestion, for every combination the majority of the predicted fragments were smaller than 500 bp (Fig 1A). In extreme cases, the Pst I- ApeK I and HinP1 I- ApeK I combinations produced a high proportion of short fragments (< 100 bp). Although the size-selection step could theoretically filter out short fragments, it is difficult to remove them completely in practice. Specifically, too many short fragments will lower the quality of library construction and subsequent sequencing. In addition, by comparing the results of three methylation-sensitive enzymes (EcoR I, HinP1 I, and ApeK I), we noticed that HinP1 I (G/CGC) could not completely digest the chicken genomic DNA in 12 h, and the ApeK I (G/CWGC) digestion products exhibited a few discrete bands (Fig 1B). By contrast, the sizes of the EcoR I (G/AATTC) digestion product were appropriate (100 bp—1000 bp) and evenly distributed without discrete bands.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Results of both in silico analysis and empirical evidence of enzyme digestions.

A) Fragment size distribution obtained by in silico digestion of the chicken genome with different double-enzyme combinations. B) Single-enzyme and double-enzyme digestion for 2 h or 12 h.

https://doi.org/10.1371/journal.pone.0179073.g001

For a more accurate estimation, we carried out direct sequencing after digestion. We presented another four parameters, including the fragment consistency index (FCI), the coefficient of variation of sequencing depth (per fragment) across three samples (CV_depth), the number of SNPs, and the distribution uniformity of SNPs across chromosomes. We prepared and sequenced 24 libraries of the eight double-enzyme digestions (three replicates for each combination, details are described in the Methods). Sequencing of all libraries produced a total of 365 million clean sequencing reads and 273 million good barcode reads, and all 24 barcode sequences were represented. The raw sequence reads were deposited in the SRA database (SRR3951559). A high FCI value represented high consistency and low level of missing data in different samples. We noticed that the fragment counts of the pooling sample (3-plex) was larger than the counts of each single sample at a high sequencing depth, indicating variance across different individuals (Table 1), which, in turn, might lead to missing data for the population. Moreover, we found that FCI was mainly determined by the types of combinations and independent of sequencing depth (FCI: 0.89, average depth: ~29× in EcoR I-Mse I; FCI: 0.64, average depth: ~36× in Pst I -Msp I; FCI: 0.58, average depth: ~16× in HinP1 I- Mse I).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Statistics of sequenced three samples from different combinations.

https://doi.org/10.1371/journal.pone.0179073.t001

The consistency of sequencing depth (per fragment) across samples is also important because it is related to genotyping accuracy. We defined the CV_depth to evaluate the performance for each of the eight combinations. The distribution of CV_depth for all fragments in each combination was shown in Fig 2. EcoR I-Mse I had the lowest mean CV_depth across three individuals (0.42±0.34 (SD)) followed shortly after by Bgl II—ApeK I (0.44±0.43 (SD)). The highest mean CV_depth across three individuals occurred in HinP1 I—ApeK I (0.82±0.51 (SD)) followed by HinP1 I—Mse I (0.77±0.51 (SD)).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Distribution of CV_depth for each combination.

Distribution of CV_depth for all fragments in each combination is displayed via boxplot. Lower and upper boundary lines of the boxes represent the 25%/75% quantile of CV_depth, and the central lines indicate the median of the data. The upper and lower whiskers represent scores outside the middle 50%. The number in each box indicate the mean of CV_depth for all fragments±standard deviation (SD).

https://doi.org/10.1371/journal.pone.0179073.g002

The number of SNPs was also critical, since too many (such as Pst I—Mse I) would increase sequencing cost, while too few (such as EcoR I—Msp I) would lower the resolution (Table 1). Another important factor for ddGBS was the distribution of SNPs per chromosome (Table 2). SNP density and coefficient of variation of SNP density (CV_{SNP density}) across different chromosomes is shown in Table 2. The SNP discovered by EcoR I-Mse I and Bgl II-ApeK I was more evenly distributed across chromosomes with a CV_{SNP density} of 0.19 than other combinations. In contrast, the highest CV_{SNP density} was found in Pst I-Mse I (0.60) followed by HinP1 I-ApeK I (0.59). Overall, the parameters among selected combinations were summarized in Table 3, and we concluded that the EcoR I- Mse I digestion was the optimal combination for the chicken ddGBS.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Distribution of SNPs discovered from three individuals across chromosomes.

https://doi.org/10.1371/journal.pone.0179073.t002

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Summary of parameters among different combinations.

https://doi.org/10.1371/journal.pone.0179073.t003

Determine the optimal level of sequencing depth

In order to obtain the optimal sequencing depth, we resampled a series of incremental subsets from total sequencing reads, and then investigated the relationship among the fragment counts, sequencing depth, and good barcode read counts. Different from whole genome sequencing, in which depth was calculated as the total length of the raw reads divided by the fixed total length of reference genome, in ddGBS the depth would be calculated as the total good barcode read counts divided by the fragment counts, which would increase with amount of sequencing until saturation. We estimated the two parameters: ORP (the minimum sequencing cost per fragment unit, Fig 3A) and SRP (the point at which there is a tangent line with zero slope for the fragment counts curve in Fig 3B) for all eight double-enzyme digestion libraries.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Relationship of ORP/SRP with good barcode reads for eight double-enzyme combinations.

A) The function of unit sequencing cost of fragments was calculated by plotting sequencing depth versus fragment counts. The ORP was defined as the minimum value of the unit sequencing cost (the minimum value of the black-dashed line). B) The sequencing reads of three individuals were sampled at five thresholds (10%, 20%, 50%, 80%, and 100%, respectively). The sequencing depth (green) was equal to the good barcode read numbers divided by the fragment counts. The SRP was the corresponding good barcode reads when the slope of the fragment counts (orange curve) reduced to zero.

https://doi.org/10.1371/journal.pone.0179073.g003

Fig 3A shows that the unit cost firstly decreased and then increased with the increasing of fragment counts. The unit cost fell to the lowest level (ORP) at approximately 10× sequencing depth in most combinations. Fig 3B shows how the fragment counts change as a function of good barcode reads. In this study, no SRPs were reached in most combinations, except EcoR I- Msp I and Bgl II- ApeK I, even though the sequencing depth was greater than 20× (Fig 3B). Theoretically, the fragment counts would be saturated when good barcode reads continued to increase, and saturation is expected to ensure consistency among different individuals. However, a typical GWAS/GS examines several hundred individuals. It is impractical to sequence all of the individuals to the saturation level (ranging in depths from 30× to 50×) for large populations (such as a family-based population with individuals of more than 100). One affordable design is to reduce the amount of sequencing appropriately and impute the missing genotypes. Taking the EcoR I- Mse I digestion as an example, at ORP, there were 2.7 million good barcode reads and approximately 270 K fragment counts for each sample, and the average sequencing depth was 10× for each fragment. Therefore, this compromise formula is not only highly precise, but also cost- effective.

SNP discovery and distribution

A total of 827 samples (824 chickens, among which three individuals were duplicated) in AIL were used to construct the ddGBS libraries. 96-plex samples were sequenced in one lane according to the ORP of EcoR I- Mse I. The raw sequence reads were deposited in the SRA database (SRR5462540, SRR5462541, SRR5462542, SRR5462543, SRR5462544, SRR5462545, SRR5462546, SRR5462547, and SRR5462548). On average, 3.44 million good barcode reads were obtained for each sample, and the average sequence depth was approximately 10×. The coefficient of variation (CV) of read counts among individuals was 0.13 (S1 Fig), indicating good consistency of library preparation. The SNPs ranged from 220–270 K among individuals prior to imputation. After strict parameter filtering in the TASSEL-BEAGLE-GBS pipeline (including imputation), we identified 291,772 SNPs ultimately (average sequencing depth was 10× with no missing data), corresponding to 1 SNP per 3.68 Kb in the chicken genome (S2 Fig and S4 Table). It is worth noting that the marker density is higher than what was reported in previous studies in chickens [16,21,22]. Among all discovered SNPs, 102,304 (accounting for 35.06% of all SNPs; the distribution is shown in S4 Table) are novel to the NCBI chicken dbSNP database (data from ftp://ftp.ncbi.nih.gov/snp/organisms/chicken_9031/VCF/ on May 4, 2016). In addition, the markers were evenly distributed without interference from GC islands and repeat regions (S2 Fig). The majority of SNPs identified were located in intergenic regions (45.19%) or intronic regions (39.55%). The exonic regions contained only 1.37% of SNPs (Table 4), comprising 51.69% missense, 3.57% nonsense, and 44.74% silent mutations.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. Number of SNPs by region.

https://doi.org/10.1371/journal.pone.0179073.t004

Genotyping accuracy evaluation

The average sequencing depth of our experiments was 10×, with an abundance of low coverage SNP sites (Fig 4A). In order to guarantee high-quality genotyping of the founders prior to imputation, we filtered raw SNP data using sequencing depth and genotyping quality, as well as minor allele frequency. Sufficient depth at each locus is essential to accurately distinguish heterozygous and homozygous sites. To assess the accuracy of genotyping, Illumina 60K chicken BeadArray microarray data and GBS results from 22 same F₀ individuals were compared. The correspondence between the two methods was evaluated at different depths ranging from 2× to 12×. When the sequencing depth reached 5×, the genotyping consistency for homozygous loci, heterozygous loci and total SNPs was 100%, 97.2%, and 99.1%, respectively (Fig 4B). The missing rates of SNPs with the 5× sequencing depth are shown in Fig 4C and 4D. After depth filtering of 5×, about 45.7% SNPs contained ≤ 50% missing genotypes (Fig 4C). The missing rates of most of the samples are between 40% and 60% (Fig 4D). We also performed two technique repeats for SNP calling for three samples, and found that reproducibility (in the case of the 5× filter condition) reached 98.5%, 98.2%, and 98.1%, respectively. Therefore, our results indicated that the genotyping results of our methods are highly reliable and accurate.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Genotyping accuracy evaluation according to different sequencing depths.

A) Distribution of coverage depths for all SNPs. B) Consistency of ddGBS genotyping results compared with 60K BeadArray microarrays using various filter conditions (sequencing depth ranging from 2× to 12×). C) Comparison of the missing rates of all 292 K SNPs on a per-site basis before and after depth filtering of 5×. D) Comparison of the missing rates of all 824 samples on a per-individual basis before and after depth filtering of 5×.

https://doi.org/10.1371/journal.pone.0179073.g004

Selection of enzyme combinations

The selection of enzymes constitutes one of the key steps in the GBS method, and is often neglected. In this study, we investigated five parameters obtained from in silico/ in vitro digestion or sequencing. We found that the result of HinP1 I enzyme digestion was not consistent with its in silico simulation. Specifically, the fragments produced by HinP1 I- Mse I and HinP1 I- ApeK I were far fewer than the results of in silico simulation. This difference could be mainly due to DNA methylation in chicken genome. Indeed, some methylation sensitive enzymes cannot digest the genome completely, which will not only cause inconsistency with the predicted results, but also interferes with the reproducibility of downstream genotype calling among different individuals [8]. Therefore, a pre-assessment experiment is necessary when markers in a new species are to be developed.

The distribution of SNPs obtained by different enzymes was tested in this study. We noticed that GGA10-20 has two-fold more fragment density than the GGA1-10 for Pst I and Msp I in silico prediction (S3 Table), which may explain the result that the fragments generated from Pst I- Msp I were not evenly distributed between macrochromosomes and microchromosomes. An uneven distribution of SNPs may have hampered the construction of evenly-distributed genetic linkage maps. However, Fábio suggested that Pst I would be suitable for chicken methylation analysis since the microchromosomes are enriched for high CpG regions [22].

Library construction process

There are two key points in library preparation that need to be addressed. First, we improved the original GBS approach described by Poland et al. in the size-selection step [23]. In previous study, 37% reads were discarded in the data processing step since the fragment size was too short (< 50 bp). Here, we removed long fragments (> 650 bp) by adding 1.3× sample volume of Agencourt^® AMPure^® XP Reagent, and removed short fragments (< 300 bp) by adding 0.8× sample volume of Agencourt^® AMPure^® XP Reagent (details are described in S1 Protocol). In our experience, magnetic bead purification was more convenient, and exhibited better consistency among different libraries compared to gel extraction.

Second, the accurate concentration of double-stranded DNA could improve the consistency of good barcode read numbers for each sample. De Donato reported that the read number per sample varied by 39% when 47 individuals were digested by Pst I [18]. Several other studies have also observed high CV in the number of different individuals (0.69 for 252-plex in chickens [16], 0.89 for 96-plex in Drosophila [39], etc.). Liao suggested that this should result from poor DNA quality, such as inaccurate quantification or contamination of DNAs with phenol/chloroform [16]. To ensure the uniformity of DNA concentration, high-molecular-weight DNA concentration was measured by Qubit2.0 prior to enzyme digestion in this study. Sequencing results showed that all of the 824 samples were well represented, and the CV of good barcode reads was 0.13 (S1 Fig), which was better than that achieved in previous studies [9,16,18,39,40]. Moreover, we noticed that the majority of missing rates of samples ranged from 40% - 60% (Fig 4D). This phenomenon may be attributed to the low CV of good barcode reads among samples.

Characteristics of sequencing depth per site

The number of SNPs declined with the increase of minimum depth used for identifying SNPs (Fig 4A), which was similar to other studies that used the GBS method [41,42]. A possible cause for the distribution was the inconsistency in the depth per fragment. The number of fragments and tags (which refers to a unique sequence from one or more good barcode reads) was counted with the Tassel parameter “-c” of 3 [33], which required a tag to be presented at least three times to be reported. Thus, a number of fragments of low depth were discarded in single sample analysis, but still retained in 3-plex pooling sample analysis, which might be the main reason for FCI < 1 (Table 1). We also noticed the divergent FCI value from the eight combinations and its dependence of combination type rather than the sequencing depth. The possible reason is that methylation (such as HinP1 I) or polymorphism disrupted a restriction site in varying extents among different samples [8]. Moreover, the difference in enzyme activity may affect the efficiency of enzyme digestion.

The accuracy of genotyping constitutes another key aspect of GBS technology. Currently, a typical GBS combines low-depth sequencing (for some inbred line of maize, the depth < 1×) and missing data imputation [43–45]. However, this strategy works better for populations with a low level of heterogeneity, such as recombinant inbred lines (RILs) rather than outcrossing populations. Genotyping errors in calling heterozygotes as homozygotes are quite common in GBS, either due to the low depth of reads or the incorrect read alignment resulting from paralogous regions. Our results showed that the 5× depth was the lowest depth for accurate SNP calls prior to the imputation step. Under this condition, the reliable SNPs of each individual were approximately 150 K, which was still higher than that in other studies [16,21]. The refined identity-by-descent (IBD) method implemented in Beagle 4.0 achieved a better performance [35].

SNP density and cost

In this study, we developed a high-density and accurate SNP genotyping method for chickens using EcoR I- Mse I. The SNP density was approximately 290 K, and some minichromosomes were not included. AIL population was commonly used for QTL fine-mapping in animal genetics [28,46–48]. Applying the SNP markers that we identified on our chicken population, we noticed that r² in the F₉ generation, which was r²_0.1 = 3.1 Kb, was substantially lower than the F₀ generation (r²_0.1 > 50 Kb) (Fig 5). This suggested that, although eight generations of recombination decrease LD levels effectively, our marker density (SNP/3.68 Kb) can still capture almost all recombination events. Thus, it can greatly benefit the fine-mapping of the QTL locations and functional genes.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. LD decay of the advanced intercross chicken lines.

A squared allelic correlation coefficient (r²) against the distance between the SNPs in the F₀ generation (HQLA was depicted as the green line, and HB was depicted as the red line) and F₉ generation (the blue line).

https://doi.org/10.1371/journal.pone.0179073.g005

The preliminary experiment was introduced to determine the optimal enzyme combination (EcoR I- Mse I) and approximate scope of the sequencing depth (ORP for large-scale population and SRP for a small number of samples). Consequently, we achieved a balance between the density of SNPs and cost. This pre-assessment method is recommended for any novel species.

Reducing cost constitutes a primary aim of all reduced-representation genome sequencing methods. Thus, our method has been optimized for cost at almost every step. However, Illumina HiSeq X Ten pair-end 150 sequencing is now much more cost-effective than single-end sequencing using NextSeq500. Moreover, adopting pair-end sequencing would provide a better chance in SNP identification than single-end sequencing in this study. Currently, ddGBS costs $30 per sample (approximately 300 K SNPs/individual), and more than 65% of the expense comes from the sequencing step in our protocol. Therefore, decreases in the cost of GBS are expected with the rapid development of sequencing technology. For example, HiSeq X Ten systems can output 800–900 Gb data/2.6–3 billion reads in a single flow cell with a cost of $1,000 for 30× of the human genome, which would be more appropriate for large-scale population sequencing (http://www.illumina.com/systems/hiseq-x-sequencing-system/system.html). The RRGS process will be quickly standardized with the declining cost of sequencing and it will, together with SNP chips, continue to be a crucial method for genomics study. In addition, the combination of RRGS and other genome-wide sampling sequencing, such as RNA-seq or Targeted Re-sequencing, could effectively promote genetics and evolutionary studies. In conclusion, we present an accurate, high-density, and cost-effective genotyping method for chickens. Our method could facilitate functional gene mapping and molecular breeding of agricultural animals, and could easily be applied to any other species.