Apple (Malus domestica) is an out-bred Rosaceae species in the Maloideae subfamily. Like most other woody tree species, its large size, cultivation area, and long juvenile period lasting up to 6 years makes it difficult to conduct genetic studies. The large chromosome number (2n = 34) with the estimated genome size is 742.3 Mb and the impossibility of producing inbred lines further complicate the matter. Despite these difficulties, there are definite advantages in studying the genetics of apple because genetically identical individuals could be obtained by vegetative propagation and grown in different environments and observed over several years (Grattapaglia and Sederoff 1994; Weeden et al. 1994; Liebhard et al. 2003a, b, c).

Genetic linkage maps allow the genomic structure to be studied, permitting the identification of quantitative trait loci (QTLs), localization of genes of interest, and providing a framework to understand the biological basis of complex traits. Construction of genetic linkage maps enables marker-assisted breeding and selection and it has proven to be especially useful for perennial tree crops such as apple, where many important traits are expressed only after years of costly field maintenance. Using marker-assisted breeding and selection, the presence of favorable alleles can be determined at the seedling stage, thereby drastically reducing the population size at the early stages of selection (Grattapaglia and Sederoff 1994; Liebhard et al. 2003c).

Several research groups have constructed high-density genetic linkage maps from both heterozygous and homozygous groups. Next-generation genotyping-by-sequencing (GBS) is now recognized as a viable single nucleotide polymorphism (SNP) validation method. GBS provides reduced representation sequencing by using restriction endonucleases to digest the genome into fragments. Using methylation-sensitive endonucleases, repetitive regions of the genome can be reduced. GBS holds great potential for plant researchers because it can be used for simultaneous marker discovery and genotyping at a reasonable cost and with a simple molecular biology workflow (Elshire et al. 2011; Spindel et al. 2013; Glaubitz et al. 2014; He et al. 2014). GBS has already been applied to the several plant species, such as barley and wheat (Poland et al. 2012), grapevine (Barba et al. 2014), rice (Spindel et al. 2013), blackcurrant (Russell et al. 2014), alfalfa (Li et al. 2014), Miscanthus sinensis (Ma et al. 2012), raspberry (Ward et al. 2013), and apple (Gardner et al. 2014). However, most of the GBS studies have focused on inbred lines partly because the commonly used bioinformatics TASSEL-GBS pipeline was designed for homozygous groups, and genotyping in heterozygous groups is more complex. To addresses this problem other GBS pipelines for heterozygous groups like raspberry (Ward et al. 2013), grape (Barba et al. 2014), and apple (Gardner et al. 2014) have been developed. Gardner et al. (2014) reported a GBS study in apple using a unique pipeline, but they mentioned that it is evident that computational and statistical tool development required to maximize analysis of GBS data. In the present study, we present diverse computational and statistical tools with the name of software and detail explanation for each step that can be used to analyze GBS data.

Here, we applied GBS methods to construct a genetic linkage map that consisted of 1,053 SNP markers over 17 linkage groups, and verified the potential application of GBS for heterozygous groups using a number of computational and statistical approaches. This novel GBS pipeline for heterozygous groups will be useful for genetic diversity analyses, phylogenetic studies, and genome wide association studies.

Plant material and DNA extraction

A segregating population from the cross between ‘Hongro’ and ‘Alps otome’ was used for the construction of the genetic linkage map. Young leaves of 93 individuals were selected from a progeny of 146 individuals developed at the Apple Research Station, National Institute of Horticultural and Herbal Science, Gunwi, Korea. Young leaves were collected and stored at −80°C until use. DNA was extracted through a modified cetyltrimethylammonium bromide method (Weeden et al. 1994) and quality was determined by agarose gel electrophoresis. DNA quantity was assessed using a SmartSpec Plus spectrophotometer (Bio-Rad, USA); each sample was diluted to 50 ng/μL.

The fruits of ‘Hongro’ and ‘Alps otome’ are highly dissimilar; ‘Hongro’ fruits are large (300 – 350 g), while ‘Alps otome’ fruits are small (25 – 50 g). Therefore, a high-quality genetic linkage map of these two cultivars will offer great potential to study the genetic control of fruit quality by QTL analysis (Kenis and Keulemans 2005; Kenis and Keulemans 2007).

Typically, heterozygous individuals are generated as a hybrid between two relatively homozygous lines. However, in apples and many other outcrossing species, this initial step is impossible. Each apple variety is highly heterozygous owing to the clonal nature of the crop and the poor performance associated with inbred material. Instead, a pseudo testcross design is typically employed in which the variety of interest is crossed to a standard variety known not to segregate for the traits being investigated (Liebhard et al. 2003c).

Genotyping-by-sequencing

For quality control of the DNA samples, we submitted the well-labeled gel image(s) to the Institute of Genomic Diversity, Cornell University. Initially, we ran 100 ng of each DNA sample on 1% agarose gels along with 500 ng of two λ HindIII size markers per gel until the standard size bands were clearly separated. Then, 300 ng of selected DNA samples (10% of samples from each plate) were digested with HindIII, a restriction enzyme that is non-methylation sensitive, and the digests were run on 1% agarose gels along with the λ HindIII size marker. Using 1 unit of restriction enzyme per sample, the DNA was digested for 2 hours at 37°C. These gel images allowed easy assessment of the quantity and quality of DNA samples (Elshire et al. 2011).

GBS was carried out at the Institute of Genomic Diversity, Cornell University. Briefly, DNA samples from the 91 F₁ individuals and the parents were digested individually with ApeKI, a methylation-sensitive restriction enzyme that recognizes a degenerate 5-bp sequence (GCWGC, where W is A or T). Then, a 95-plex GBS sequencing library was prepared by ligating the digested DNA to unique barcoded adapters and common adapters. PCR amplifications were performed according to a standard PCR protocol (Elshire et al. 2011). Sequencing was performed on an Illumina HiSeq 2000 platform (San Diego, CA). The parental DNAs were sequenced twice to increase the accuracy of the called SNPs.

Preprocessing

After sequencing, the raw reads were de-multiplexed according to the barcode sequences and trimmed using a Python script (Fig. 1). This script splits the raw Illumina FASTQ file into 93 separate FASTQ files based on the barcode sequences associated with each sample and filters out reads that contain any ambiguous bases in the barcodes. The reads that contained only the common adapter were also trimmed using the Cutadapt software (Martin 2011).

The de-multiplexed reads were trimmed using the Solexa QA package v.1.13 (Cox et al. 2010). It is common for the quality of bases to decrease at the ends of Illumina reads, so we trimmed the ends when the Phred quality score dropped below Q = 20 (or 0.05 probability of error). Additionally, we trimmed all 5′ and 3′ stretches of ambiguous N nucleotides. Poor quality sequence reads and reads shorter than 25 bases were discarded.

Alignment and detection of SNPs and InDels

We used the Burrows-Wheeler Aligner (BWA; 0.6.1-r104) program (Li and Durbin 2009) to align the clean reads to the apple (Malus × domestica) Whole Genome v1.0p as the reference genome (Velasco et al. 2010). The BWA default values were used for mapping, except for seed length (−l) = 30, maximum differences in the seed (−k) = 1, number of threads (−t) = 32, mismatch penalty (−M) = 6, gap open penalty (−O) = 15, and gap extension penalty (−E) = 8. Mapped reads were extracted for further analyses from the resulting BAM file using SAMtools v.0.1.16 (Li et al. 2009). The high mapping quality ensured reliable (unique) mapping of the reads, which is important for variant calling. Using the SAMtools varFilter command, SNPs were called only for variable positions with a minimal mapping quality (−Q) of 30. The minimum and maximum of read depths were set to 3 and 100, respectively. An in-house script that considered biallelic loci was used to select the significant sites in the called SNPs positions (Kim et al. 2014). Depending on the proportion of SNP/InDel reads to mapped reads, the variant types were classified into three categories: homozygous SNP/InDel for more than have 90% of same SNP at loci, heterozygous SNP/InDel for more than 40% and less than 60%, and other for those that remained.

SNP filtering and segregation analysis

SNP calls were filtered for quality by restricting the marker set to biallelic SNPs, and requiring genotype calls at each SNP to have a depth of three reads in each sample. Missing SNP allele data in individuals were imputed using a Perl script. The script was used to identify and genotype valid and high-quality SNPs in the 91 F₁ mapping individuals along with the two parental genotypes using a sliding window approach. Haplotype phasing involved sliding a window along the chromosome, estimating haplotype phases within each window and piecing fragments together over the whole chromosome. The window size was set according to a major allele frequency of 0.7.

The marker code was determined using JoinMap 4.0 (Van Ooijen 2006) and population type CP (cross pollinators). CP indicates a cross between two heterozygous diploid parents, with linkage phases originally unknown (Van Ooijen 2006). Five segregation types of CP populations (lmxll, nnxnp, hkxhk, efxeg, abxcd) have been described (Van Ooijen 2006), but only three of the segregation types (lmxll, nnxnp, hkxhk) were genotyped in this study. Segregation type ‘nnxnp’ describes markers with first parent being homologous and second parent being heterozygous, ‘lmxll’ describes markers with first parent being heterozygous and second parent being homologous, and ‘hkxhk’ describes markers with both parents being heterozygous. Then, valid loci for genetic mapping were filtered using the following criteria. First, the expected segregation ratio for ‘lmxll’ and ‘nnxnp’ is 1:1, and including ‘hkxhk’ is 1:2:1. Only markers whose segregation ratios did not significantly deviate from the expected ratios were retained. The pattern of allelic segregation was tested using a chi-squared (χ²) test to determine deviations from expected Mendelian segregation ratios, and SNPs with significant segregation distortion (χ²P-value < 0.01) were excluded from further analysis. Second, the sequencing depth of each locus of each plant was checked, and genotypes with low depth (> 3 for homozygous SNP, > 5 for heterozygous SNP) were set as missing, because they could be caused by sequencing or other errors. Third, any locus from nnxnp and lmxll with more than 10 missing data was filtered, because too much missing data could influence mapping results. However, the criteria for hkxhk was more than 20 missing data, because it has relatively smaller number than other segregation types.

Linkage map construction

To increase mapping efficiency, pairs or groups of loci with identical genotypes (i.e., complete linkage) were identified and a single marker was chosen to represent the group. Linkage analysis was performed with JoinMap 4.0 using the parameters set for the regression mapping algorithm. Linkage groups were estimated by applying independence LOD (Logarithm of Odds) threshold ranges from 2 to 13, and constructed with a linkage LOD of at least 8.0, a recombination fraction of 0.35, a map LOD value of 1, a goodness-of-fit jump threshold of 3, and a ripple value of 1. The recombination frequencies were converted into map distances (cM) using the Kosambi mapping function (Kosambi 1943).

Genotyping-by-sequencing

A total of 223 million sequence reads were obtained from the mapping population (Table 1). To ensure sufficient read depths at potential SNP loci and to increase the accuracy of the called SNPs, we performed deeper sequencing of the parental lines, obtaining 4,332,022 and 5,013,265 raw reads for ‘Hongro’ and ‘Alps otome’ respectively. We obtained reliable trimmed data by applying the Python script with strict criteria (probability value = 0.05, Phred quality score = 20, and minimum read length bp = 25 bp). From the standard TASSEL-GBS pipeline, we obtained 29 million sequence reads, which was insufficient to build an apple genetic linkage map (Table 1). However, from KNU-GBS pipeline, we obtained nearly 170 million sequence reads with an average of 1,827,937 reads per sample and used this trimmed data for further analysis (Supplementary Table S1, Table 1).

Alignment of trimmed GBS reads to the reference genome

A total of 64.08% of the trimmed reads were mapped to the reference genome, the average depth of the mapped regions was 5.37, and 2.25% of the reference genome was covered (Supplementary Table S1). The trimmed reads from the parental samples covered a relatively high percentage of the reference genome (7.7% and 7.79% for ‘Hongro’ and ‘Alps otome’ respectively), most likely because of the deeper sequencing of the parental lines. We also verified that the trimmed reads were well distributed on the reference genome by plotting the average depth of the mapped regions to the positions (Mb) of the mapped reads on the reference genome (Fig. 2).

Recombination pattern analysis in the segregation type

We verified the reliability of our filtering criteria by analyzing the segregation type ‘nnxnp’ from chromosome 13 (Fig. 3). Before imputation, the segregation types consisted dominantly of ‘nn’ and a recombination pattern could not be determined. However, after conducting the strict imputation, the natural ratio of the ‘nn’ and ‘np’ types with clear recombination regions was observed. Although some regions still showed an unnatural recombination pattern, we presumed that these regions could potentially be recombination hotspot regions in chromosome 13.

Genetic linkage map construction

Total numbers of 547,131 SNPs were detected using the SAMtools varFilter command and 28,761 SNP were divided into three types of segregation type: 13,912 from nnxnp, 5,264 from lmxll, and 9,585 from hkxhk. Then we performed the chi-squared (χ²) test with significance level ≤ 0.01, total 9,843 SNPs (5,757 from nnxnp, 2,519 from lmxll, and 1,567 from hkxhk) were selected. As it mentioned in materials and methods, after applying different cutoff for missing data total number of 7,274 SNPs (4,612 from nnxnp, 1,775 from lmxll and 887 from hkxhk) were selected. After removing the redundant SNP markers, the final set of 2,590 SNPs remained for linkage map construction. Of these, 836 were present in the parental marker code ‘nnxnp’, 867 were in the marker code ‘lmxll’, and the remaining 887 were in the marker code ‘hkxhk’ (Table 2).

Based on the filtering criteria from JoinMap 4.0, described in the material and methods section, we constructed a genetic linkage map from the ‘Hongro’ × ‘Alps otome’ cross that consisted of 1,053 SNP markers over 17 linkage groups (Fig. 4). The map encompassed 1350.1 cM, with linkage groups ranging from 56.3 cM (LG14) to 97.2 cM (LG3). The numbers of SNP markers mapped to each linkage group varied from 34 in LG4 to 123 in LG2, with a mean of 61.9 SNPs per linkage group. The average marker density was 1.28 cM per marker and maximum gap size ranged from 5 cM in LG12 to 21.8 cM in LG17 (Table 3).

GBS data analysis and mapping

We obtained nearly 170 million sequence reads with an average of 1,827,937 reads per sample using KNU-GBS pipeline. This number is much higher than the average of 973,896 reads per sample reported previously from an apple population (Gardner et al. 2014). The difference between the two apple GBS analyses may be influenced by the characteristics of the mapping populations. This is because; they used the ‘Golden Delicious’ variety that was used to sequence the apple reference genome, as the maternal line for the mapping population (Velasco et al. 2010).

F₁ population mapping with GBS provides many advantages such as simple, less laborious, and highly reproducible (Elshire et al. 2011). Furthermore, repetitive regions that are prevalent in apple genome, can be avoided using methylation-sensitive REs. For this reason, it is evident fact that choice of the restriction enzyme is also one of the key points for GBS. We used the ApeKI restriction enzyme, a 5-base cutter that cuts frequently, and the Gardner group used the double restriction enzymes, HindIII-HF (6-base cutter) and MspI (4-base cutter). Based on their choice of restriction enzymes, the Gardner group may have obtained more fragments per sample. However, it is more complex to conduct the GBS with double restriction enzymes than one enzyme. In our research, we verified that even using one ApeKI restriction enzyme in apple, it would be enough to make dense genetic linkage map.

Inversions of SNP order

Several inversions of SNP order from their expected physical locations were observed from the GBS analysis. One reason for this observed order inversion is the characteristics of the reference genome. For example, the apple genome contains a high frequency of repeated regions; repetitive elements correspond to 500.7 Mb (67%) and 98% of the unassembled genome sequences are repetitive (Velasco et al. 2010). These high percentages of repetitive genomic regions can lead to incorrect anchoring of sequences to the apple reference genome. For example, Gardner et al. (2014) also reported that many SNPs on LG9 were assigned to other chromosomes such as chromosome13 and 17 according to the apple reference genome.

Genetic linkage map construction

Although the number of SNPs in our final set was higher than in the previous study (Gardner et al. 2014) (2,590 and 2,436 respectively), after filtering, the number of SNPs that remained in the genetic linkage map was smaller than in the previous study (1,053 and 1,994 respectively). This smaller number is a result of the different LOD thresholds that were used for map construction: a linkage LOD in our map is 8.0 compared with a threshold of 6.0 in the previous study (Gardner et al. 2014) to increase the accuracy of the map. Even more stringent LOD score, 9.0, was used in pear (Wu et al. 2014). For this reason, the average marker density of our group is also smaller than the previous study (1.28 cM and 0.68 cM respectively). Since our study was focusing on SNPs, InDels what we found here would be an interesting extension of this work.

In this study, we describe the construction of a genetic map for M. domestica consisting of 17 linkage groups based on 1,053 SNP markers spanning 1,350.1 cM using a new developed GBS pipeline, not only because existing TASSEL-GBS pipeline is mainly designed for homozygous group. We demonstrated that GBS methods could be applied successfully to heterozygous species for the construction of a high-density linkage map. We also show how the utility of the GBS data can be maximized by applying the computational and statistical tools in our newly developed GBS pipeline. This pipeline and analyses with constructed genetic linkage map reported here could become a valuable tool for QTL analysis of agronomically important traits, marker-assisted breeding programs, and comparative and structural genome analyses.

This work was supported financially by the Next-Generation BioGreen21 Program (PJ01311503), Rural Development Administration, Republic of Korea.