A high-density, multi-parental SNP genetic map on apple validates a new mapping approach for outcrossing species

Abstract

Quantitative trait loci (QTL) mapping approaches rely on the correct ordering of molecular markers along the chromosomes, which can be obtained from genetic linkage maps or a reference genome sequence. For apple (Malus domestica Borkh), the genome sequence v1 and v2 could not meet this need; therefore, a novel approach was devised to develop a dense genetic linkage map, providing the most reliable marker-loci order for the highest possible number of markers. The approach was based on four strategies: (i) the use of multiple full-sib families, (ii) the reduction of missing information through the use of HaploBlocks and alternative calling procedures for single-nucleotide polymorphism (SNP) markers, (iii) the construction of a single backcross-type data set including all families, and (iv) a two-step map generation procedure based on the sequential inclusion of markers. The map comprises 15?417 SNP markers, clustered in 3?K HaploBlock markers spanning 1?267?cM, with an average distance between adjacent markers of 0.37?cM and a maximum distance of 3.29?cM. Moreover, chromosome 5 was oriented according to its homoeologous chromosome 10. This map was useful to improve the apple genome sequence, design the Axiom Apple 480?K SNP array and perform multifamily-based QTL studies. Its collinearity with the genome sequences v1 and v3 are reported. To our knowledge, this is the shortest published SNP map in apple, while including the largest number of markers, families and individuals. This result validates our methodology, proving its value for the construction of integrated linkage maps for any outbreeding species.

Introduction

Genetic linkage maps play a major role in clarifying the genetic control of important traits and the development of DNA-based diagnostic tools for marker-assisted breeding. They are supposed to reflect the order of genes and molecular markers as they occur on the chromosomes and are critical resources for: (i) the identification of gene location on chromosomes via quantitative trait loci (QTL) discovery studies,^1–3 (ii) the building of reference genome sequences through anchoring, ordering and orienting of contigs and scaffolds,⁴ and (iii) the cloning of genes through map-based approaches.^5–7 Most of the economically important traits in plant breeding, such as yield and product quality, are quantitative and controlled by multiple genes. Therefore, identifying the genomic location of such genes is a high priority for selecting new improved crop varieties.^8,9 Remarkable advances have been achieved in understanding the functional complexity underpinning quantitative traits. A number of QTL with strong effects on phenotypic variation have been discovered, genetically positioned, validated and, in various cases, successfully exploited in marker-assisted breeding.^9–11

In outbreeding species, conventional QTL discovery approaches rely on the availability of genetic linkage maps and segregating bi-parental full-sib (FS) families. However, a single FS family is unlikely to segregate for all QTL, thus providing only partial information. Currently, QTL mapping is shifting toward the simultaneous analysis of more complex pedigreed FS-families, derived by multiple direct parents and founders.^12–18 This approach increases the probability of detecting QTL and capturing allelic variation while it improves the characterization of QTL performance in different genetic backgrounds.^12,19–21

The EU-funded FruitBreedomics project¹⁰ was aimed, among other objectives, to clarify the genetic determination of a series of fruit quality traits in apple through a multifamily QTL mapping approach using molecular markers from a 20?K Infinium SNP array.²² This raised the need for a reference genetic linkage map allowing adequate integration of SNP marker data across wide germplasm. The accurateness of marker order is crucial to remove sources of spurious double recombinants and to narrow the intervals where QTL are located. When a high-quality consensus map or reference genome sequence is available, they can be used for the correct ordering of markers.

At the onset of this work, various genetic linkage maps were available for apple with most based on a single FS family^2,23–32 and some based on a few FS families.^33–35 Furthermore, a draft apple reference genome sequence was available,³⁶ which has been used for developing whole-genome genotyping (WGG) assays^22,36–38 for producing high-density SNP linkage maps on segregating FS families.^28,34,35 However, all of the array and Genotyping By Sequencing (GBS) derived genetic linkage maps highlighted discrepancies in marker positions to the reference genome for ~14 (ref. 28) to 22% (ref. 32) of the markers. The generation of a highly reliable integrated map was prompted by the inconsistencies among these maps and by their low proportion of common markers.

Genetic linkage maps are created through the study of co-segregation patterns of markers and genes in segregating families. In outbreeding species, usually both parents can contribute segregation information and the generation of three different linkage maps is allowed: two parent specific maps and one integrated bi-parental map. A relevant issue in the construction of integrated genetic maps on bi-allelic markers, such as SNPs, is the high proportion of non-informative data, which are due to three main causes. First, missing values are inevitably high as most of the SNPs segregate in only one parent, thus being homozygous and not informative for the second parent. Second, markers segregating in both parents (abռab?aa, ab, bb) yield only 50% of informative data since the alleles of the ab progeny genotypes cannot be unequivocally traced back to the donor parent. The reduction in information is even worse when a null allele is present in both parental genotypes (a-ռa-): since their progeny is called by the presence (a- and aa) and absence (--), 75% of the genotypes (a- and aa) cannot be unequivocally called and will be uninformative. Third, most markers usually do not segregate in each family. Therefore, the total amount of missing information goes well beyond 50% for any SNP marker. Uncertainty in marker order may also arise from standard approaches for map integration when merging the two parental maps of a FS family into a single bi-parental map through automated procedures.^39,40 This process raises ambiguity in the appropriateness of marker order due to incompleteness in segregation information. Accordingly, linkage map integration across multiple bi-parental maps further increases ambiguities due to rise in missing data.

The main purpose of the present work was to produce a highly reliable and high-density integrated multi-parent genetic linkage map for apple (Malus domestica Borkh) to be used as a reference genetic map and as support in improving the apple genome assembly. To obtain the most reliable order for the highest possible number of markers, a novel mapping procedure was adopted by combining the following four main strategies: (i) using 21 segregating FS-families genotyped with the recent 20?K Infinium SNP array;²² (ii) reducing the proportion of non-informative data through an ad hoc SNP filtering and calling method and by the use of the HaploBlock (HB) bins formed by tightly linked markers; (iii) using a backcross (BC) design on single-parent data, rather than a Cross-Pollinator (CP) design on bi-parental ones, to facilitate the integration of parental data (full details are explained in the methods section); and (iv) using a two-step mapping procedure where an Initial Framework Map (IFM) of only the most informative markers, provides a reliable starting point for adding the remaining less informative.

Materials and methods

Plant material

This study included 1586 progeny plants from 21 FS families.²² They were obtained from 26 parents and originated from six different breeding programs from five European countries (Table 1). Eighteen of them were part of the previous European project HiDRAS.⁴¹ Although most of the families comprised ~50 individuals, 7 of them significantly differed in size (Table 1), with 12_J being the smallest (23 individuals) and DLO.12 the largest (219 individuals). Part of these populations have also been used in studies on the development of multiplexes of SSR markers,⁴² validation of the pedigree-based analyses approach for QTL mapping on multiple pedigreed families¹³ and QTL discovery for horticultural traits.¹⁸

DNA isolation

For each individual, young, preferably unfolded leaves were sampled and freeze dried. Genomic DNA was extracted according to Schouten et al.²⁷ The DNA was further purified using an RNAse step and quantified using 0.8% agarose gels and a dilution series of an external reference Lambda DNA (Invitrogen, Carlsbad, CA, USA).

SNP-genotyping

The 21 FS families and their parents were genotyped with the 20?K Infinium SNP array at the Fondazione Edmund Mach according to published procedures.^28,37 SNP calling was performed using GenomeStudio software (Illumina Inc., San Diego, CA, USA; http://www.illumina.com), with a GenCall threshold of 0.15, and ASSIsT,⁴³ a filtering and calling pipeline that accounts for the presence of null-alleles and signal intensity differences among AB-genotypes, thus increasing the number of usable SNP markers and providing some fully informative bi-parental markers with segregation type abռac.

SNP markers origin and focal point (FP) design

The 20?K Infinium SNP array consists of three different SNP sources: the recently designed FruitBreedomics (FB) markers²² and subsets of the RosBREED SNPs and GD_snp markers (jointly referred to as IRSC—International Rosaceae SNP Consortium-markers) present on the previous 8?K Infinium SNP array.^22,37 The GD_snp markers are based on polymorphisms within Golden Delicious sequence data, which were previously validated using SNPlex technology.³⁶ The new FB markers were designed in clusters of up to 11 SNPs located within a region of at most 10?kbp, defined as Focal Point (FP) and distributed along the genome at distances of ~1?cM.²² The 8?K Infinium SNP array also followed an FP design for the IRSC markers; however, here each FP stretched up to 100?kbp.^22,37

Construction of bi-parental single-family linkage maps and SNP validation

Before generating the integrated multi-parent linkage map, integrated bi-parental SNP linkage maps were created for each of the 21 FS-families to validate filtered SNPs and to verify the tight linkage between SNP markers coming from the same FP. For the construction of genetic linkage maps, JoinMap v4.1 (ref. 44) was used, applying the multipoint maximum likelihood mapping algorithm approach for cross pollinators^40,45 and the Haldane mapping function using pre-set default settings. Markers were removed from the data set of an individual FS-family when they showed a severely distorted segregation (P<0.01) and nearby markers segregating for the same parent did not show such a distortion, or when the rare genotype class occurred in less than 5% of the progeny. The GenomeStudio cluster plots were examined for the following: (i) markers mapping far (>10?cM) from any other marker; (ii) markers showing high nearest neighbor stress (NN stress) values (>2?cM) according to JoinMap output, and (iii) genotype calls involved in double recombinant single-points (singletons) as reported by JoinMap. When necessary and feasible, calls or the parental origin of markers were adjusted. Markers that remained problematic were excluded. Markers with identical genotypic scores (identicals), which are automatically set aside by JoinMap, were added back to the resulting linkage maps.

SNP markers from the same probe that mapped on different LGs in distinct families were classified as multi-locus SNPs, and they were considered as distinct markers specifying the mapping LG in their name. Also, the ‘$’ symbol is added in front of their initial name to easily distinguish them from single-locus SNPs. Moreover, whenever SNP markers from the same FP were mapping to distinct genetic regions, the FP was split into two or more distinct SNP clusters, or in cases of an individual SNP, this was moved out of the FP. Finally, the assessed sets of co-segregating SNP markers belonging to the same FP were defined as HBs.

SNP assignment to HBs

Validated SNPs, as mapped in at least one of the 21 FS families (Supplementary Table S1), were grouped into HBs according to the FP design on the 20?K array^22,37 or to their coordinates on the reference genome sequence v2 (in case of IRSC markers)²² while accounting for their co-segregation patterns in the individual FS families. HBs comprising only FB markers (FB-HBs) spanned at most 10?kbp, those including only IRSC markers (IRSC-HBs) covered at most 100?kbp in size, and a window of 20?kbp was assigned to HBs that comprised both FB and IRSC markers (FB+IRSC-HBs). Those SNPs that did not fall in any physical distance range allowed by the FP design were set aside the HBs and kept as individual SNP markers. Within each HB, SNP markers were ordered according to the coordinates of the targeted sequence polymorphism on the above mentioned genome sequence.

The HB marker and the BC strategy

The creation of HBs of co-segregating markers allowed a bin-mapping strategy where the segregation information of adjacent SNPs was aggregated and condensed into a single, virtual HB marker. The aggregation of co-segregating markers within the same HB increases the genotype score robustness consequently to information redundancy, and marker informativeness when combining markers with different segregation types. This is the case when a marker segregating in a single parent (for example, abռaa) is combined with a bi-allelic marker heterozygous in both parental plants (for example, abռab) or with a single parent marker of the other sex (for example, aaռab), leading to the generation of a fully informative marker record (corresponding to a segregation type abռac, or abռcd). In view of our mapping effort, this strategy was implemented in the ad hoc developed software Haploblock Aggregator (HapAg—http://www.wageningenur.nl/en/show/HaploblockAggregator.htm) and applied to our data (Supplementary File S1). For each FS family, HapAg aggregated the segregation information of the SNP markers belonging to the same HB by using the information on linkage group and the linkage phase of the individual markers (Figure 1a), while considering the meiotic events occurring in the two parental plants separately. Thereto, HapAg splits the parental allelic contribution of every individual of each family into two distinct sub-data sets including either maternal or paternal recombination events (Figure 1b). Eventually, maternal and paternal data sets of all the progenies from all FS-families were merged to generate a single BC-type data set (Figure 1c), having twice the number of individuals as the original CP populations (a more detailed description of the methodological steps performed by HapAg is available in the software manual (http://www.wageningenur.nl/en/show/HaploblockAggregator.htm)). The BC segregation type allows the correct phasing of the markers segregating in different families, leading to integrating the genotypic data prior to map construction and to the production of a unique integrated genetic map rather than a map resulting from the a posteriori integration of the linkage maps obtained from FS families.

When inconsistencies between SNP scoring within the same HB are present (for example, due to recombination or gene conversion within the HB or calling issues), warning messages are reported, and the aggregated data score is set to missing (--) by HapAg. The warnings were carefully examined to identify the origin of the issue and, when possible, to solve it before producing the final data set. Specifically, when conflicts occurred in more than five individuals per HB across all families, these were always inspected, while less consideration was given to a lower number of warnings as compatible to the expected number of true double recombination. Given the overall estimated genetic to physical size ratio of 548?kbp/cM, based on the results from Velasco and colleagues,³⁶ the probability of a recombination event occurring within a single HB is expected to be 1.8ձ0^-4 (FB-HBs), 1.8ձ0^-3 (IRSC-HBs) and 3.5ձ0^-4 (FB+IRSC-HBs). These probabilities are so low that we decided to neglect such recombination, at least during the construction of the map, and thus, make missing the HB marker call for the potentially recombinant individual.

The final BC-type data set, formed by HB markers and individual SNPs not included in HBs (for simplicity, the generic term ‘markers’ will be used to refer to both), was used to construct the integrated Genetic Linkage Map (iGLMap).

iGLMap construction

Linkage maps were produced using JoinMap v4.1 with the same algorithm and parameter settings used when mapping the single FS-families, but now following a BC design. Due to the sensitivity of mapping algorithms to missing data,⁴⁰ a two-step mapping procedure was adopted. First, a reliable IFM was determined by including only markers that had genotypic data for at least 25% of the individuals (from 800 to ~3?200 available meiosis). Subsequently, the map was completed by adding less informative markers with up to 90% missing data (from 320 to 800 available meiosis), while using the obtained IFM marker order as the Start Order in JoinMap. Markers with more than 90% missing values (327) were not considered. At both steps, for each linkage map, double recombinant single points (singletons) identified by the JoinMap Genotype probabilities function were visually examined through a graphical genotyping approach⁴⁶ whereby data were displayed in map order as color-coded genotypes in Microsoft Excel using the conditional cell formatting feature. Issues causing singletons were investigated. For each LG, the best map was defined as the one with the lowest number of singletons.

The iGLMap for apple genome improvement

The produced IFM was used to examine and improve the reference apple genome sequence from v2, as used for the design of the 20?K Infinium SNP array,²² to v3 (https://www.rosaceae.org/species/malus/malus_x_domestica/genome_v3.0.a1). Specifically, this map helped to assess potential issues in the anchoring of contigs within scaffolds and provided a new assignment of scaffolds to chromosomes.

To confirm the collinearity between v3 of the reference genome and the final iGLMap, physical coordinates of all genetically mapped SNP markers were plotted as a function of genetic distances on the iGLMap. For comparison, collinearity between the iGLMap and the most widely used apple genome sequence v1 was also evaluated. MareyMaps⁴⁷ for each chromosome were produced using R (R Development Team Core, Vienna, Austria, Europe).

Results

Construction of 21 bi-parental linkage maps

The 21 FS families (Table 1) were genotyped using the 20?K Infinium SNP array.²² Overall, 15?697 SNP markers (87%) passed the SNP calling and filtering pipeline ASSIsT⁴³ and were genetically mapped using a final set of 1?586 individuals. Single-parent map integration produced high-quality bi-parental maps for each family and LG (Table 2; Supplementary File S2), with six exceptions. LG6 and LG16 in family I_W and LG7 in family I_CC show large homozygous regions which led either to the absence of segregating bi-parental bridge markers (thus only single parent maps could be generated) (Supplementary Figures S1a and b), or to the lack of segregating markers in one parent (LG7 I_CC). Furthermore, LG17 in family I_CC has no integrated map because the few bi-parental markers are not sufficiently spaced to allow orientation of the two single parent maps. This family is also peculiar since its paternal LG17 shows a highly distorted segregating region, coincident with the self-incompatibility S locus (Supplementary Figure S1c), thus suggesting the two parents to have a common self-incompatibility allele. Finally, the LGs 13 in families 12_J and 12_K have a genetic length of 1.5–1.6x the population mean (Table 2). Such relative large sizes are usually an indicator for data issues. This was true for 12_J, where the family size showed to be too small to come for a meaningful genetic map for this LG. For 12_K however we suspect a biological reason as its distorted segregations and double recombination pattern may be explained by the presence of natural selection at distant loci of the same chromosome (data not shown).

The total length of the 21 maps ranges from 1123 (GaPi) to 1?551?cM (12_K) with an average length of 1?305?cM (±113.3) across families (Table 2a). The average distance between SNPs ranges from 0.16?cM in JoPr to 0.25 in 12_K. On average across families, 6?848 (±619) SNPs were mapped ranging from 5?570 (12_N) to 8?454 (JoPr) (Table 2b).

Focal points validation and HBs identification

The order of markers on the 21 genetic maps was thoroughly checked, and correspondence to the apple’s physical map v2 was assessed. Overall, the SNP marker order was coherent. Nonetheless, the following three main issues were identified that were common to all families and, overall, affected 22% of the markers: (i) inconsistencies in LG assignment for 17% of the markers; (ii) regions of inversion (max ~2.5?Mbp) involving 3% of the markers; and (iii) misplaced regions of markers within the same pseudo-chromosome for 2% of the markers. Furthermore, a close inspection of the bi-parental genetic maps highlighted the presence of 111 multi-locus SNPs (0.7% of the markers), mapping in 2 or even 3 (on 4 occasions) different LGs across distinct families (Supplementary Table S2); 29 of these fall in known homoeologous regions originating from an ancient apple genome duplication.³⁶ The 111 multi-locus SNP probes resulted in additional 115 alternative mapped SNP loci, establishing a final total number of 15?812 validated SNP-loci, each segregating in at least one family (Table 3, Supplementary Table S1).

Therefore, the expected tight linkage between SNP markers of the same FP could not always be confirmed because 91 FPs included SNPs mapping to distinct genetic regions. Such discrepancies led to the creation of independent SNPs sub-clusters, identified as distinct HBs or individual SNPs. At this stage, 921 individual SNP markers and 2?837 HBs were considered (Table 4).

HB marker genotype definition and BC data set

The software HapAg was designed in the framework of the current mapping effort to aggregate the segregation information of the SNPs belonging to each HB into a single HB marker. During the aggregation process, HapAg reported 2?848 inconsistencies in SNPs scoring within 748 HBs (26%). More than half of them (51%) resulted from the erroneous assignment of 55 SNPs to a FP and, therefore, to a HB, due to inadequate genome sequence information (unidentified during the mapping effort on single families). Those SNPs were removed from the HBs and used as individual SNP markers (Supplementary Table S1). The remaining 49% of conflicts involved 25% of the HBs and were due to calling errors, recombination within HBs, or gene conversion.

The majority of these HBs (96.3%) presented no more than 5 warnings, which were considered having a minor impact. The remaining 3.7% reported between 6 and 16 conflicts. These inconsistencies were addressed by setting the HB score of the conflicting individual as missing. In fact, it would have been unfeasible to examine all cluster plots discriminating between calling errors, true recombination events, or gene conversion in a reasonable amount of time.

The final integrated data set consisted of a single BC-type population of 3?172 individuals, including the genotypic information of 2837 HB makers and 976 individual SNPs. The overall proportion of missing information was massively reduced from 78% of the initial set of 15?812 SNP markers to 54% of the final (HB+individual SNPs) integrated data set (Figure 2), thus retaining the complete genetic information of the larger SNP data set.

Construction of the iGLMap

An IFM was constructed using markers that had data for at least 25% of the 3?172 individuals, which was true for 2631 HB and 344 individual SNP markers. The genetic length of its 17 LGs ranged from 64 (LG1) to 113?cM (LG15), for a total IFM length of 1?279?cM (Table 5). The IFM was subsequently completed by adding markers with genotypic data for at least 10% of the individuals, to produce the iGLMap. Over the whole mapping process, an improvement in map quality was obtained through the close inspection of 1?320 singletons by graphical genotyping (Supplementary File S3), whereby 67% of the singletons showed to result from 386 misplaced markers. These were moved to nearby positions. The adjusted marker order was then used as the Start Order for map re-estimation verifying that singletons’ issues were indeed solved, and no new double recombinants were generated by the shift. Next, 13% of the singletons came from 13 HB markers that could not be adequately placed along the map, because they showed conflicting genotypes with adjacent markers at any position. Since the underlying cause could not be identified, those HB-markers (Supplementary Table S1) were removed. Finally, 20% of the initial 1?320 singletons remained. These are likely caused by genotyping errors, gene conversion events or true double recombinations. As the examination of each individual case would have been too time intensive, singletons’ scores were set to missing, following Bassil et al.⁴⁸

The final iGLMap (Figure 3, Supplementary Figure S4) utilized 15?417 SNPs (Table 3) and consisted of 3?473 markers as follows: 2?797 HBs (Table 4), composed of a total of 14?741 SNPs, and 676 individual SNP markers (Table 5). The total genetic length is 1?267?cM, with LG1 being the shortest linkage group (63?cM) and LG15 being the longest (112?cM). The maximum distance between adjacent markers of 3.29?cM is observed on LG6 (Table 5). An estimate of marker position robustness is provided by JoinMap Plausible positions output table where, for each marker, the probability of its assigned and alternative positions are estimated across 1000 iterations. Graphical results for each LG are presented in Supplementary File S4. The order of markers appeared highly reliable, as all 1000 iterations resulted in a single position for most markers. Only a few regions show plausible alternative marker orders. Generally, such regions are characterized by very high marker density. When only two markers are involved, the distance between such markers is usually less than 0.01?cM. When three or more markers are involved, they usually span a region of less than 0.05?cM (with the only exception of four regions in total, one on each of the linkage groups LG1, LG9, LG10 and LG15, which spanned from 0.07 up to 0.1?cM). Thus, the overall results show a robust and accurate marker order along the iGLMap with just a few small regions of uncertainty where the resolution was not enough to distinguish the genetic order of tightly linked markers.

iGLMap collinearity with the apple genome

The 21 bi-parental genetic maps already led to the identification of many regions where the genetic and physical maps deviate from collinearity. The final iGLMap allowed for visualization of inconsistencies with version v1 and v3 of the apple genome³⁶ (Figure 4). The most used version is v1, which consists of a Primary Assembly and various alternative scaffold sequences (Figure 4a) representing homologous and possibly also some homoeologous sequences. Genetic data suggests the presence of some inverted regions that can be recognized by short lines at a right angle to the main curve of each graph in Figure 4a, which goes from the bottom-left to the top-right. Two examples are approximately at 20?cM on LG10 and at 35?cM on LG11. The short sequences of dots clearly aligned above or below the main curve represents regions where the position on the genome sequence may be shifted (that is, approximately at 80?cM on LG10, at 75?cM on LG11 and at 60?cM on LG14). Some regions of the iGLMap have their sequence counterparts on non-corresponding pseudo-chromosomes. For example, the initial part of LG1 and the 22–34?cM region of LG10 include markers belonging to other nine pseudo-chromosomes according to the genome sequence; also, the initial region of LG4 (0–10?cM) consists almost entirely of markers from the pseudo-chromosome 9. Alternative scaffold sequences may also combine information from multiple pseudo-chromosomes. For example, the alternative scaffold on pseudo-chromosome 9 seems to contain, for the 20?Mb region, a small fragment of its homoeologous region on LG17.

To support the development of a new 487?K Axiom array,³⁸ version 3 of the apple reference genome was developed by using, among others, genetic information from the IFM. Adjustments were made only when entire contigs could be moved. Though overall collinearity has greatly improved, 2% of the genetically mapped markers is not represented anymore, and inconsistency is still observed for 6% of the markers. As shown in Figure 4b, collinearity still varies considerably as follows: some LGs need further attention for just some narrow regions (for example, LG1 and LG17), while others have several regions to be ameliorated (for example, LG13 and LG16). At least six LGs show lines parallel to the main curve, suggesting the presence of scaffold sequence segments that have been misaligned (LG7, LG9, LG13, LG14, LG15 and LG16).

Discussion

The iGLMap was developed through an innovative mapping approach whose power and major novelties rely on the use of HBs and a single BC-type data set. The HB approach provided highly informative common markers across families, facilitating integration of the maps. The use of the BC strategy applied to an outcrossing species allowed the integration of the genotypic data prior to map construction, avoiding the multiple efforts of parental map integration within and across FS families. The advantages of the HB strategy and BC design were strengthened by the use of an improved SNP filtering tool (ASSIsT⁴³), which increased the number of informative markers; and by a multi-step map construction, which facilitated the inclusion of markers with relatively high missing value content in the final high-density linkage map.

Use of multiple FS families

The production of multifamily genetic maps can lead to an increase in marker density and genome coverage, overcoming local loss of genetic resolution and thus having the potential for a more accurate order of markers.⁴⁹ To our knowledge, the iGLMap features the highest number of families ever used for a crop plant’s genetic mapping, followed by 13 families for wheat.⁴⁹ Indeed, the merging of data from 21 FS families resulted in an integrated map with an average of 2.3 times the number of markers of its preceding single family maps and a genetic length highly consistent with that of the preceding maps (284?cM shorter than the longest and 144?cM longer than the shortest one). The enhanced map resolution is indicated by the consistent reduction in the maximum distance between adjacent markers as follows: 3.29?cM on the LG6 of the iGLMap, compared to values for individual families ranging from 5.94 (LG7 of the FuGa family) to 30.14?cM (LG13 of the 12_K family). The use of multiple families thus allowed for overcoming the limitations that occurred in single families including issues related to large regions of homozygosity and to extremely skewed segregation ratios, as those reported here for families I_W and I_CC, and solved in the final iGLMap.

Use of improved SNP calling and HBs

The essence of our strategy is the reduction of missing information, which usually hampers accurate marker ordering. Therefore, we aimed to include as many common SNPs as possible by exploiting information from null-alleles, and to increase information content of bi-parental abռab markers also by exploiting signal intensity differences.^43,50 This contrasts with standard SNP calling procedures, which usually do not recognize null-alleles and discard such markers while classifying them as ‘non-Mendelian’.

The use of null-alleles is expected to improve mapping efficiency not only by increasing the number of common markers, but also by increasing information content of single SNP markers: this allows the use of tri-allelic markers of type a0ռb0 that provide full segregation information for both parents.^43,50 Accounting for differences in signal intensity due to an additional SNP at the probe site also gives rise to fully informative markers, because two variants may occur from one or both of the marker alleles converting bi-allelic abռab markers into tri-allelic abռab’ ones.^43,50

Our data set included 1?244 SNP markers with null-alleles that segregated in at least one FS-family, resulting in additional data for a total of 1?974 marker/family combinations. Of these, 327 markers resulted in 494 (25%) combinations that showed full segregation information (abռac), 514 markers out of 915 (46%) combinations that segregated as dominant abռab’ markers for which a null-allele was present in both parents. Additionally, 328 initial abռab markers became fully informative abռac markers for 471 combinations by exploiting signal intensity differences resulting from an additional SNP at the probe site. The accounting for null-alleles and for differences in signal intensity thus increased the amount of segregation information for HB-markers and individual SNP-markers. The impact of this method on the final data set will vary with the study population and will increase at decreasing numbers of families and HB size.

The proposed HB strategy may show some similarities with bin mapping, although it presents substantial differences that can lead to a more accurate marker positioning. Usually, bin mapping is based on the identification of an interval (bin) along a linkage group where, given a small set of individuals, no recombination events are observed for any of them.⁵¹ This approach has been very successful to obtain approximate genetic information on marker position along a chromosome using a small number of individuals.^52–54 However, the smallest unit of resolution in these genetic maps was given by the bin, whose length usually spanned from an average size of 7.8?cM up to a maximum length of 33?cM.^52,54 Conversely, our HB strategy defined bins based on (very short) physical distances. The use of multiple SNPs within a narrow physical bin was anticipated in the design of the 8 and 20?K Infinium SNP arrays through the identification of tightly clustered SNP markers within a small physical window.^22,37

Indeed, the HB strategy was very effective as it led to an increase in marker robustness through the exploitation of redundant genetic information of adjacent markers, and to an increase in information content, especially when markers were of different segregation type. In fact, in pseudo-testcross analysis,³⁹ bi-allelic markers heterozygous in both parents lead to 50% non-informative data as it is not possible to unequivocally establish the origin of the alleles in the heterozygous progenies. On the contrary, the proposed HB strategy is able to exploit these data to produce a robust and fully informative HB marker when such bi-parental markers are combined with markers of other segregation type(s) within the same HB (see family 2 in Figure 1). Consequently, whereas all single bi-allelic SNP markers were informative for less than 50% of the progenies across multiple families (Figure 3), 55% of the aggregated HB-markers are informative for more than 50% of the individuals, and 8.6% of the HB-markers presented 90 to 100% of the information, thus providing highly informative bridge markers common to most FS-families. The use of HB-markers thus favored map integration across families and allowed high accuracy in marker position along the chromosome.

In this study, we did not further examine within-HB SNPs order, and we disregarded a priori any recombination inside HBs by making the relative marker scores missing. However, the within-HB recombinations may still be useful in the future to validate local genome assemblies when needed. We estimate that ~1?940 within HB recombinations will be present in the whole data set. This estimation is based on a genome-wide relationship between physical and genetic distances of 586?Kb/cM (the estimated genome size of 742?Mb [ref. 36] divided by the length of the iGLMap), the physical length covered by the HBs (using the number of mapped HBs and their maximum allowed size for each of the three different HB-type), the total number of individuals (3?172), the availability of informative data (46%), and assuming absence of crossover interference.

This estimated number is certainly inflated, as crossover interference does occur. Furthermore, the distance between the two most apart and informative SNPs within a HB is less than the allowed maximum HB size. Also, not all these recombinations will be noticed as having occurred within the HB, because, in many parents, HBs lack the two or more informative SNPs needed to observe recombination. Consequently, the number of expected observable recombination events might not be very different from the 1?396 conflicts observed among SNP calls within HBs. Notably, if such estimations are correct, they also imply the absence of major calling issues in the final data set.

The HB approach makes use of all the available markers, including identical markers sensu JoinMap. These are markers that have exactly identical genotype calls across all progenies of a single mapping family, including an identical linkage phase. In the analysis of single families JoinMap removes such identicals from the data set by default, thus reducing memory requirements and computation time. However, we re-entered them because they are needed for the construction of HB markers. The availability of full genetic information for each HB, completed by the use of the identicals, allows the adequate mapping of these HBs. Furthermore, in the case of conflicting data within a HB, the presence of identical markers may be helpful in tracing the inconsistency.

BC mapping strategy

The most commonly used linkage analysis approach for outbreeding FS families is based on the two-way pseudo-testcross strategy,³⁹ which results in the production of two single-parent maps to be then integrated using the available bi-parental bridging markers. Despite the advances in map integration methods^40,55–58 the process is still based on computing the average of marker distances over the two parental maps, which is affected by the different segregation types and by the number and informativeness of bridging bi-parental markers. The consequent result is a loss of parental-specific features, low accuracy in marker order⁴ and inconsistencies between individual maps. Similarly, the most popular approaches used in the generation of multi-population consensus maps^59–61 are based on the integration of marker distances over populations. Therefore, map integration is not a straightforward process, especially because usually not all the loci are common in all populations, and additional factors, such as the local reordering and marker misplacement, further affect resolution and accuracy of the consensus map.

On the contrary, the proposed BC design, leads to the generation of a single data set where genotypic data have been merged prior to the generation of integrated genetic maps. This strategy reduces the integration process to a single step where segregation data, recombination events and marker order are directly related and accessible for close inspection across all germplasm through visual approaches, such as graphical genotyping.⁴⁶ The use of this strategy to obtain high-quality linkage maps is demonstrated by the short length of the iGLMap and the robustness of its marker order (see below). Moreover, the use of a condensed, fully informative data set optimizes computational performance.

A multi-step mapping procedure and data curation through visual inspection

Multi-step procedures are very often used for the construction of high-density consensus maps as they allow for marker order optimization and genotyping error correction.⁶² In standard mapping approaches, bi-parental segregating framework markers can be used to produce back-bone maps, which provide a fixed order on which the final consensus map is built by adding the remaining markers.⁶³Bin-mapping methods also follow a two-phase mapping process based on the initial construction of high-confidence framework maps to which new markers are subsequently added.^25,51,64

In our case, the choice to pursue a two-step mapping procedure was further motivated by the sensitivity of the JoinMap mapping algorithms to the presence of missing data.⁴⁰ In fact, JoinMap was not able to perform the grouping when analyzing all 3473 markers simultaneously. For this reason, we first decided to produce a robust IFM including only the most informative markers.

Nevertheless, in high-density marker data sets, the marker order proposed by JoinMap demonstrated sensitivity to the presence of single problematic data points⁴⁸ and/or structured missing data as present in our final data set due to markers segregating in some families but not in others. This may lead to inadequate marker orders that can cause inflation of map size, a high number of singletons and larger regions of double recombination. In Bassil et al.,⁴⁸ as in our current study, the use of a graphical genotyping tool was crucial to perform adequate data curation through identifying problematic data points and alternative marker orders for improving the map quality and reducing singletons and double recombination.

Evaluation of map quality

High quality and robustness of the iGLMap is inferred from its genetic length when compared to other published linkage maps of apple. Since the release of the genome sequence of apple³⁶ and the consequent extensive development of whole genome genotyping SNP arrays,^22,37 several high-density genetic maps of apple were constructed. Most of them are integrated bi-parental maps of one family, where the number of individuals ranges from 118 to 297, comprising different type of markers (for example, SNP+SSR for a total of 1,016 markers;³⁴ SNP+SSR for a total of 2?579 markers;²⁸ SNP+SSR for a total of 601 markers;⁶⁵ 3?441 SNP (ref. 66)) and GBS based linkage maps comprising 3?967 SNP (ref. 29), 1?053 SNP (ref. 31), and 1?918 and 2?818 SNP markers (ref. 32). Furthermore, a consensus map of apple was produced by merging a SNP genetic map with maps from earlier studies on eight FS families³⁴ as follows: two based on SSR markers^24,25 and six based on SSR and SNP markers.³⁶ This consensus map consisted of 2?875 markers with a total genetic length of 1?991?cM. Finally, a multifamily single parent map was produced for the ‘Honeycrisp’ apple from three small progenies (318 individuals in total) having 1?091 SNPs (ref. 35). All of these genetic maps have a total size ranging from 1?282–1?991?cM. ^28,34 The iGLMap is shorter than all these despite having the highest chance of false recombination due to data issues; and therefore an inflated genetic length due to being based on the largest number of markers, the largest number of FS families and the largest total number of individuals. This result once more proves the quality of the iGLMap and thereby the validity of our approach.

Further support for the quality of the iGLMap comes from the plots of the alternative plausible positions calculated by JoinMap. In the iGLMap, spots with uncertain marker order usually spanned less than 0.05?cM, and only occasionally spanned up to 0.1?cM. Such limited uncertainty was feasible thanks to the reduced level of missing values, the high level of data curation, and the high number of individuals that were included in our study, which theoretically allow a maximum resolution of ~0.02?cM when two HB markers segregate in all individuals.

The length of the iGLMap (1?267?cM) is shorter than both the direct and weighted average length of the 21 single family maps (1?305?cM), because additional data curation had been performed after the construction of the individual maps, during the data aggregation and the mapping of the obtained HB data, as described in the result section.

Proposal for a re-orientation of the LG5 linkage map

In apple and pear, the convention for orienting linkage groups and pseudo-molecules is based on the genetic linkage map by Maliepaard et al.² for ‘Prima’בFiesta’. Maliepaard et al.² were the first to report large stretches of homoeologous sequences as detected by multi-locus targeting RFLP clones. However, they did not use these first indications on homoeology for linkage group orientation. To date, those initial observations have been confirmed, and knowledge about duplication patterns have considerably increased.³⁶ In apple, the orientation of the two chromosomes of the homoeologous pairs 5 and 10 and 13 and 16, where homoeology involves almost the entire chromosome, showed opposite orientation sensu Maliepaard et al.² The orientation of LG13 has been modified to match that of LG16 by Liebhard et al.⁶⁷ In the present work, we modified the orientation of the LG5 linkage map to match that of LG10 to facilitate future chromosome and genome comparison studies between homoeologous chromosomes within a species and synteny studies among species. We re-orientated LG5 and not LG10 because LG5 has been reported in less genetic studies on QTL and candidate genes identification in apple (Supplementary File S5). Due to the limited number of cases where putative QTLs and genes have been reported on LG5, we expect the benefits for future studies to offset the initial inconvenience generated by comparing our iGLMap with those of previous papers. In pear, the naming and orientation of linkage groups was based on apple⁶⁸ having made use of the high synteny between these two crops.^68,69 In the few QTL reports on pear, LG5 has been mentioned more frequently than LG10 (Supplementary File S5).

The HB strategy, genotyping by sequencing and single families

The HB strategy contributed considerably to the quality of the IGLMap. This strategy is applicable to a wide range of genetic data and population types. Indeed, in this study, HB-sizing was essentially based on information from a draft genome sequence: Haploblock markers were created by the aggregation of segregation information from SNP that co-localized within a narrow physical window. The approach is also applicable with GBS derived data. In case of, Restriction site Associated DNA (RAD) sequencing,⁷⁰ the various SNP from the same read may be used for the construction of a HB marker. The HB strategy may also be used on single families as it may considerably reduce computation time, especially in case of a large sized family and very high marker densities.

Conclusion

The reliable marker order, high coverage and resolution of the iGLMap makes it a valid reference map for QTL mapping studies and the evaluation of genome assemblies. Therefore, the iGLMap can contribute to the enhancement of marker-assisted breeding approaches aimed at improving apple quality and productivity. The usefulness of the iGLMap has already been demonstrated in at least the following two occasions: i) its successful use in a recent pedigree-based analyses study identifying QTLs and possibly the underlying candidate genes for controlling chilling and heat requirements¹⁸ and ii) the improvement of the reference genome sequence of apple (https://www.rosaceae.org/species/malus/malus_x_domestica/genome_v3.0.a1).

Finally, the methodology presented here may be valuable for the construction of accurate high-density bi- and multi-parental integrated genetic linkage maps for any outbreeding species.

Competing interests

The authors declare no conflict of interest.

Acknowledgements

We thank Yolanda Noordijk for the isolation of DNA from all samples at Wageningen-UR and Elisa Banchi for her work on the genotyping of these samples with the 20?K Infinium SNP array at the Fondazione Edmund Mach. This work has been co-funded by the EU seventh Framework Programme by the FruitBreedomics project N°. 265582: Integrated Approach for increasing breeding efficiency in fruit tree crops (www.FruitBreedomics.com). The views expressed in this work are the sole responsibility of the authors and do not necessary reflect the views of the European Commission.

References

Footnotes