High-resolution genetic maps of Lotus japonicus and L. burttii based on re-sequencing of recombinant inbred lines.

Recombinant inbred lines (RILs) derived from bi-parental populations are stable genetic resources, which are widely used for constructing genetic linkage maps. These genetic maps are essential for QTL mapping and can aid contig and scaffold anchoring in the final stages of genome assembly. In this study, two Lotus sp. RIL populations, Lotus japonicus MG-20 × Gifu and Gifu × L. burttii, were characterized by Illumina re-sequencing. Genotyping of 187 MG-20 × Gifu RILs at 87,140 marker positions and 96 Gifu × L. burttii RILs at 357,973 marker positions allowed us to accurately identify 1,929 recombination breakpoints in the MG-20 × Gifu RILs and 1,044 breakpoints in the Gifu × L. burttii population. The resulting high-density genetic maps now facilitate high-accuracy QTL mapping, identification of reference genome mis-assemblies, and characterization of structural variants.

1. Introduction

Legumes are important in both agricultural and natural ecosystems because of their capacity for symbiotic nitrogen fixation, and many are rich sources of protein., The legume Lotus japonicus has been extensively used as a model for deciphering the molecular genetics governing the symbiotic interaction with rhizobia., Comprehensive genetic and genomic resources have been developed in L. japonicus, including bacterial artificial chromosome libraries,, a LORE1 mutant population,, a TILLING population, ESTs,, genetic and physical maps, and recombinant inbred line (RIL) populations.

So far, L. japonicus genetic maps have been based on low-resolution simple sequence repeat (SSR) genotyping of F2 (∼900 markers) and RIL populations (∼96 markers)., These maps have proved useful for anchoring large contigs to pseudomolecules and for QTL mapping,, but both approaches have been limited by the low marker densities of the SSR-based maps. The resolution of a genetic map depends on the number of recombination events in the mapping population and on the number of markers used to genotype each individual. The recombination breakpoint positions in the current L. japonicus maps are determined with limited accuracy due to the low marker density. With advances in sequencing technology, it has now become possible to genotype all single nucleotide polymorphism (SNP) positions in RIL populations using whole-genome re-sequencing. SNP genotyping using various methods has been successfully used in other plant species. In soybean and Medicago sativa, genetic maps were generated by genotyping ∼3,000 SNPs,, in chickpea ∼2–4,000 SNPs were used, whereas 16,000 single feature polymorphisms were used to generate a genetic map with 815 recombination breakpoints in Arabidopsis. In Brassica rapa more than 1 million SNPs were discovered in 150 RILs using whole-genome re-sequencing, and a genetic map with 2,305 recombination events was constructed. Similarly, a high-density genetic map with 5,074 breakpoints was generated in rice using low-coverage re-sequencing of 150 RILs.

The choice of crossing partners and the population size is important for generation of high-quality genetic maps. The F2 populations developed earlier using Gifu as the common parent in crosses with Funakura or L. filicaulis, suffered drawbacks including low levels of polymorphism, low viability, severe segregation distortion and recombination suppression., Although the MG-20 × Gifu cross has advantages in terms of good viability with little segregation distortion and a reasonable marker density,,, it has the disadvantage that a large reciprocal translocation between the top of Gifu chromosome 1 and the bottom of MG-20 chromosome 2 causes severe suppression of recombination in this genomic region.

In order to overcome the problems with suppression of recombination in the MG-20 × Gifu cross and to increase genetic variation and marker density, a more divergent accession originating from West Pakistan, Lotus burttii,, was used for bi-parental crossing with Gifu. F2 populations derived from this cross showed good viability and did not suffer from suppression of recombination at the top of chromosome 1, and a genetic map based on genotyping 97 SSR markers in 146 RILs has been constructed.

Here, we have genotyped 187 MG-20 × Gifu and 96 Gifu × L. burttii RILs using whole-genome re-sequencing to accurately identify recombination breakpoints and construct high-resolution genetic maps for both crosses. These maps now facilitate identification of reference genome mis-assemblies, assignment of unanchored contigs onto pseudomolecules, and mapping of genomic translocations with high accuracy.

2. Materials and methods

2.1. Sequencing and read alignment

DNA from 187 MG-20 × Gifu (F9) and 96 RILs of Gifu × L. burttii (F8) RILs developed from F2 populations by single seed descent, was extracted using the CTAB protocol. Sequencing libraries were then constructed using the Illumina Nextera DNA library preparation kit (FC-121-1031) with a dual-index adapter system according to the manufacturer’s instructions. The libraries were sequenced on an Illumina HiSeq 2000 sequencer to generate 2×93 bp paired-end reads from ∼500 bp insert libraries. The reads were mapped to the MG-20 L. japonicus genome v.3.0 (http://www.kazusa.or.jp/lotus/, 17 June 2016, date last accessed) using Burrows–Wheeler Aligner mem v. 0.7 with default parameters .5a.

2.2. Genetic map construction

Candidate polymorphic positions were available from genome sequences, and these were validated by verifying that the positions were polymorphic according to the RIL genotyping results. The candidate SNP markers were genotyped using the mpileup function of SAMtools software v0.1.19. The command used to determine genotypes was samtools mpileup –uD –b –f –I | samtools bcftools view –cg - > result.vcf’, where –u is to generate uncompressed binary variant call format (BCF) output, -D is to output per sample depth. Because the reads were mapped to an MG-20 reference sequence, the Gifu × L. burttii RIL genotype calls were re-oriented with reference to Gifu, such that the Gifu, rather than the MG-20, genotype was considered the reference allele.

To generate genetic maps, we used a modified version of a previously described approach. Because of the higher sequencing depth used in our study, we simplified the sliding window approach for calling genotype blocks, using overlapping windows with a single base offset spanning 20 markers upstream and downstream of each marker position. Genotype counts were calculated for each window, and the consensus genotype was called as the genotype with the highest abundance. Recombination breakpoints were then called at positions where the consensus genotype switched (Fig. 2B). To ensure better accuracy, the genotype block calls were compared with the original genotyping data at each marker position for each RIL and a score was calculated. The score was +1 if the genotype block call and marker genotypes matched, 0 if there was a mismatch, and +0.5 if one of the genotypes was heterozygous. This marker score was used to identify inaccurate breakpoint positions, where the breakpoint had been miscalled by a few marker positions relative to its true position. Inaccurate breakpoint positions were then shifted to optimize the marker score. For the MG-20 × Gifu data, the marker score was additionally used for filtering out markers that showed poor concordance with the genotype block calls. Markers were removed if the cumulative score at each marker position was <150. Finally, breakpoints that occurred at the same position in the majority of the RILs were removed manually because these very likely represented artifacts. In order to construct a genetic segment map for each of the RIL populations, the region between any two consecutive breakpoints in all RILs for each population was defined as one segment.

Figure 2.

Genotype block calling. Chromosome 1 is shown for the MG-20 × Gifu cross. Genotype counts in 41-marker windows is shown. (A) Unfiltered genotype cells and (B) Final genotype blocks.

2.3. Identification of translocations and assembly errors based on the Gifu × L. burttii genetic map

For the MG-20 × Gifu RILs, the segment map did not require any filtering, but for Gifu × L. burttii, a large number of aberrant segments were apparent (Fig. 3B). Although the complex segment pattern presented was difficult to resolve accurately using automated approaches, we identified aberrant segments with high confidence by manual curation of the segment map, where the dominant genotype pattern for each genomic region was readily apparent and the aberrant segments stood out because of large deviations from the consensus genotype pattern. Next, the aberrant segments that were <10 kb were discarded in order to focus on the largest and most well-supported segments. Each of the remaining aberrant segments was then matched with all other segments using a custom perl script to determine the genetic linkage. If ≥90 out of the 96 RIL genotypes were identical between two segments, then the aberrant segment was placed next to the matching segment, identifying it is misplaced in the original map. Unplaced aberrant segments were not included in the genetic map.

Figure 3.

Recombination fraction and LOD scores for all segment pairs. Recombination fractions (rf) are shown above the diagonal, and the LOD scores below. Low recombination fractions and high LOD scores indicate strong segment linkage. (A) MG-20 × Gifu RILs, (B–D) Gifu × L. burttii RILs. (B) Including aberrant segments. (C) After replacing aberrant segments. (D) After replacing aberrant segments and the chromosome 1–2 translocation.

The displaced segments were then categorized as translocations or assembly errors. The RIL genotypes of the displaced segment in the Gifu × L. burttii map were compared with the MG-20 × Gifu RIL genotypes in the same genomic region by manual inspection. If the region in MG-20 × Gifu showed markers with discordant genotypes compared with neighboring regions, then the segment in Gifu × L. burttii was classified as an assembly error. In contrast, if the genotypes of the markers in the region matched the flanking markers, the segment was categorized as a translocation. The genetic maps were visualized by calculating pairwise logarithm of the odds (LOD) scores and recombination fractions using the plot.rf function of the R/qtl package.

3. Results and discussion

3.1. RIL re-sequencing and genotyping

Paired-end sequencing reads from 187 MG-20 × Gifu and 96 Gifu × L. burttii RILs were aligned to the L. japonicus MG-20 reference genome v. 3.0 using bwa. The resulting average coverage was 6.53 for the MG20 × Gifu RILs and 8.03 for the Gifu × L. burttii RILs (Fig. 1 and Supplementary Table S1). Next, 87,140 positions polymorphic between Gifu and MG-20, and 357,973 positions polymorphic between Gifu and L. burttii were genotyped in the two sets of RILs on chromosomes 1–6 (Supplementary Files S1–S2). With ∼230 Mbp sequence anchored on these pseudomolecules, this yields average marker distances of 2.6 and 0.64 kbp in the two maps, allowing high-resolution mapping of recombination breakpoints. In addition, 64,011 positions polymorphic between MG-20 and Gifu and 203,540 positions polymorphic between Gifu and L. burttii were genotyped on 162 Mbp of unanchored contigs (chr0) (Supplementary Files S3–S4) in L. japonicus v 3.0. This genotype information can now be used to position chr0 contigs onto the L. japonicus v. 3.0 pseudomolecules.

Figure 1.

Average RIL read coverage for (A) MG-20 × Gifu and (B) Gifu × Lotus burttii.

3.2. Identification of RIL genotype blocks and construction of genetic maps

Because recombination events are relatively infrequent, we expected to identify large genotype blocks based on the genotyping data. These blocks were apparent, but they were frequently interrupted by aberrant genotype calls (Fig. 2A). This noise could be due to genotyping errors or to misplaced markers, and we used a sliding window approach to reduce the noise by calling consensus genotypes within 41-marker sliding windows (Fig. 2B). In addition, we eliminated small blocks of markers, which showed very poor consistency in genotype calls with those of their neighbors (Supplementary Files S5–S6).

Following identification of genotype blocks for all RILs, we combined the results into genetic maps, where each segment constituted the region between any two consecutive breakpoints in the RIL population. To visualize the resulting maps, we then plotted the pairwise recombination fractions and LOD scores for all segments (Fig. 3). The MG-20 × Gifu map contained 1,744 segments and appeared very consistent with previous results,, with the strongest linkage along the diagonal and at the site of the Gifu – MG-20 chromosome 1–2 translocation (Fig. 3A, Supplementary Table S3). In contrast, with 1,837 detected breakpoints, the Gifu × L. burttii map contained many more segments than expected for 96 RILs, and large numbers of aberrant linkage signals were apparent as numerous inconsistent horizontal and vertical lines in the LOD score and recombination fraction plot (Fig. 3B).

3.3. Identification of assembly errors and correction of the Gifu × L. burttii map

Because all reads were aligned to the MG-20 reference genome, the inconsistencies in the Gifu × L. burttii map could be caused either by assembly errors in the MG-20 reference sequence being detected because of the fourfold larger marker density in the Gifu × L. burttii cross or by the occurrence of translocations specific to MG-20. To distinguish between these possibilities, we returned to the unfiltered MG-20 × Gifu genotyping results (Supplementary File S1) and investigated the genotype consistencies in the aberrant regions identified in the Gifu × L. burttii map. In 179 out of 347 cases, we found markers within the same intervals in the MG-20 × Gifu map, which had been filtered away because of genotype inconsistency with neighboring markers, identifying these intervals as misplaced in the MG-20 reference assembly. This genetic information has potential for improving the quality of the reference genome assembly because an average marker distance of only 0.6 kbp in the Gifu × L. burttii map would reliably detect misplaced contigs with sizes down to a few kilobasepair.

We then proceeded to replace the aberrant segments in the Gifu × L. burttii map based on their genotype patterns, and we were able to unambiguously assign 201 out of 347 regions, resulting in a consistent map comprising 1,044 segments (Fig. 3C, Table 1) Finally, the reciprocally translocated regions at the top of chromosome 1 and bottom of chromosome 2 were replaced to align the segment map with the genetic data (Fig. 3D, Supplementary Table S4, Supplementary File S4).

Table 1.

Summary of replaced segments in the Gifu × L. burttii segment map

Chr	chr1	chr2	chr3	chr4	chr5	chr6	Total
chr1	24	9	10	8	6	1	58
chr2	12	8	6	6	2	2	36
chr3	12	6	12	3	3	2	38
chr4	6	9	1	8	3	2	29
chr5	8	1	2	1	2	0	14
chr6	2	13	3	2	3	3	26
Total	64	46	34	28	19	10	201

The segments were moved from the chromosome indicated in the header row to the chromosome indicated in the leftmost column.

3.4. Segregation distortion and breakpoint distribution

After establishing the genetic maps for both crosses, we characterized segregation distortion patterns. Although the MG-20 × Gifu population was well balanced, we observed more severe biases in the Gifu × L. burttii RILs. The most striking signal was found on chromosome 2 with a maxium at ∼25 Mbp. Here, the L. burttii allele frequency reached 96% (Fig. 4), indicating very strong selection against Gifu alleles resulting in elimination of 50% of the possible progeny due to this incompatibility alone.

Figure 4.

Segregation distortion. The Gifu genotype fractions for MG-20 × Gifu and Gifu × Lotus burttii are plotted for each segment in the genetic maps. The positions are relative to the MG-20 v. 3.0 reference sequence.

Both crosses showed very similar averages of 10.3 breakpoints per chromosome (Fig. 5, Supplementary Table S2), and the genetic length of each chromosome corresponded to the previous genetic linkage map (Supplementary Table S2). In addition, the breakpoints and overall genotype patterns detected in the MG-20 × Gifu RILs corresponded well with previous genotype data generated using SSR markers, with the exception of seven lines, RI-011, RI-068, RI-080, RI-097, RI-149, RI-152 and RI-153, where the line identities may have been interchanged.

Figure 5.

Recombination breakpoint counts. (A and B) Histograms of breakpoint counts by chromosome. (A) MG-20 × Gifu RILs. (B) Gifu × Lotus burttii RILs. (C and D) Histogram of breakpoint counts averaged across chromosomes 1–6. (C) MG-20 × Gifu RILs. (D) Gifu × L. burttii RILs. Vertical lines indicate the average number of breakpoints.

In the MG-20 × Gifu cross, extensive regions displaying suppression of recombination were located on chromosomes 1 and 6, whereas the Gifu × L. burttii cross showed suppression of recombination on chromosomes 5 and 6 (Fig. 6). Both crosses showed a reduced recombination frequency at the top of chromosome 2, which corresponds to a region rich in ribosomal RNA genes, and near the centromere on chromosome 4.

Figure 6.

Recombination breakpoint distribution. Vertical lines represent the positions of recombination breakpoints. The black dot on the grey bar indicates the centromere position for each chromosome. The horizontal bars indicate regions of recombination suppression. All breakpoint positions are shown relative to the MG-20 v.3.0 reference sequence.

3.5. Mapping of MG-20-specific structural variants

Because suppression of recombination can be a consequence of structural rearrangements and because a major reciprocal chromosomal translocation has occurred between chromosomes 1 and 2 in MG-20 with respect to Gifu and L. burttii,, we returned to the aberrant segments in the Gifu × L. burttii map in order to precisely characterize the structural variation between MG-20 and Gifu/L. burttii. Only two aberrant segments were found in regions, which contained consistent marker genotypes in the MG-20 × Gifu population (Supplementary Fig. S1). These segments corresponded to the regions displaying suppressed recombination in the MG-20 × Gifu cross, the top of chromosome 1 and a region around 20 Mbp on chromosome 6, indicating that they likely represented structural variants (Fig. 6, Supplementary Fig. S1).

Examining first the top of chromosome 1 and the chromosome 1–2 reciprocal translocation, we could accurately map the translocation breakpoints to between 1,730 and 1,806 kbp on chromosome 1 (Fig. 7A) and between 32,017 kbp and 32,645 kbp on chromosome 2 (Fig. 7C) using the Gifu × L. burttii map. A number of putative structural rearrangements were also detected within this translocated region based on the Gifu × L. burttii data, but these were all classified as assembly errors and were mostly caused by errors in scaffolding, as the aberrant segments corresponded to entire contigs created by short read assembly.

Figure 7.

Mapping structural variants. (A) Translocation breakpoint at the top of MG-20 chromosome 1. (B) Translocated segment from chromosome 1 placed at the correct location at the bottom of chromosome 2. (C) Translocation breakpoint at the bottom of MG-20 chromosome 2. (D) Translocated segment from chromosome 2 placed at the correct location at the top of chromosome 1. (E) Translocation breakpoints on chromosome 6. (F) The translocated segment placed at the correct position on Gifu chromosome 6. Number ranges indicate possible breakpoint intervals. Numbers above arrows indicate segment endpoints.

For the region on chromosome 6, genetic data from the Gifu × L. burttii map indicated translocation of a 186 kb region from ∼21.4 Mbp in MG-20 to ∼17.8 Mbp in Gifu (Fig. 7E and F). Although there was no indication of a larger inversion, it appeared that this smaller translocation could be causing the suppression of recombination in the region around 20 Mbp (Fig. 6), which was also observed in the existing genetic map of an MG-20 × Gifu F2 population, where more than 10 markers mapped to 46.6 cM, corresponding to the same region.

In conclusion, the two high-density genetic maps are complementary and together provide adequate coverage of the L. japonicus genome for QTL mapping and genome assembly, and the very high marker density generated by whole-genome re-sequencing of the Gifu × L. burttii RILs is especially useful for identifying and correcting even minor reference genome mis-assemblies.

4. Availability

The RILs are available from LegumeBase (http://www.legumebase.brc.miyazaki-u.ac.jp/). The genetic maps and all marker genotypes are provided as supplemental files and tables. The raw sequencing data are available from the DDBJ Sequence Read Archive with accession numbers DRA004729, DRA002730 and DRA004731.