The black grain barley Hatiexi No.1 was crossed with the yellow grain variety Zhe5819. The resulting F1 plants were self-crossed to obtain F2. Grain color of F1 and F2 were examined in the field before harvested. The F2 population of Hatiexi No.1 × Zhe5819 consists of 551 black grain lines and 172 yellow grain lines. For mapping the gene controlling grain color, homologous yellow individuals were selected from F2 population of Hatiexi No.1 × Zhe5819.

Young leaves of the two parents (Hatiexi No.1 and Zhe5819) and F2 individuals were collected for DNA extraction. Total genomic DNA was prepared from leaf tissues using CTAB method (Murray and Thompson, 1980). DNA concentration and quality were estimated using a Nanodrop 2000 UV-vis spectrophotometer machine and by electrophoresis through 0.8% agrose gels.

Fifty plants with black grain and fifty plants with yellow grain were selected randomly from the F2 generation as two pools for SLAF-seq-BSA. The black pool and yellow pool were constructed by mixing an equal amount of DNA from 50 black individuals and 50 yellow individuals, respectively. The parents and two pools were used for SLAF library construction and sequencing as described previously (Sun et al., 2013; Hu et al., 2016). A pre-design SLAF experiment was designed to determine conditions and appropriate restriction enzymes for digestion that optimize SLAF yield and maximize SLAF-seq efficiency. The SLAF library was conducted in accordance using the pre-designed scheme. Genomic DNA was digested with RsaI (New England Biolabs, NEB). After that, a single-nucleotide A overhang were added to the digested fragments with Klenow Fragment (3′→ 5′ exo–) (NEB) and dATP at 37°C, and then the Duplex Tag-labeled Sequencing adapters (PAGE purified, Life Technologies) were ligated to the A-tailed fragments with T4 DNA ligase. PCR reaction was performed using diluted restriction-ligation DNA samples, dNTP, Q5 High-Fidelity DNA Polymerase and PCR primers: AATGATACGGCGACCACCGA and CAAGCAGAAGACGGCATACG (PAGE purified, Life Technologies). The PCR productions were purified using Agencourt AMPure XP beads (Beckman Coulter, High Wycombe, United Kingdom) and then pooled. The pooled sample was separated by electrophoresis using 2% agarose gel. Fragments with 364–394 bp (with indexes and adaptors) in size were excised, and then purified using QIAquick Gel Extraction Kit (QIAGEN). The gel-purified product was sequenced on the Illumina HiSeq 2500 system (Illumina, Inc; San Diego, CA, United States) according to the manufacturer’s recommendations.

After sequencing, all reads were aligned to barley reference genome released by The International Barley Sequencing Consortium in 2012 (IBSC 2012¹) using BWA software (Li and Durbin, 2009). Sequences with over 90% identity were grouped in one SLAF locus. Specific fragments were considered as SLAF tags and polymorphic SLAFs were selected due to their polymorphism between two parents. Based on physical position of SLAF tags, SNP calling was performed by local realignment and mutation detection using GATK software². We excluded SNPs which supported less than four reads in the two pools and showed no polymorphism between the parents because they may be false positives due to genomic repeat sequence, sequencing or alignment errors. Then SNPs showed multiple allele loci and monomorphism between the black and yellow pools were removed. Finally, SNPs with one genotype derived from Hatiexi No.1 and the other from Zhe5819 were identified as polymorphic markers, and were selected for association analysis.

Association mapping was conducted to identify candidate regions for black lemma and pericarp using both SNP_index (Abe et al., 2012; Takagi et al., 2013) and Euclidean distance (ED) methods (Hill et al., 2013).

SNP_index association analysis, a recently published method, is used to calculate genotype frequency differences between two bulks that were satisfied by Δ(SNP_index). The closer marker is associated with phenotype while the closer Δ(SNP_index) is associated with 1. M stands for Hatiexi No.1, P stands for Zhe5819, aa denotes the genotype from Hatiexi No.1 in the black pool, and ab denotes the genotype from the yellow pool. The Δ(SNP_index) was calculated as follows:

Mab and Pab were the depth of yellow pool from black and yellow grain parents, respectively; and Maa and Paa indicated the depth of black pool from black and yellow grain parents, respectively.

The allelic frequency was calculated by Euclidean distance followed by Loess regression analysis which identifies regions in which QTL lies and generates a list of putative regions in the linked genomic segment.

Euclidean distance association analysis is a type of method that calculates Euclidean distance and is satisfied by ED according to Hill et al. (2013) and Geng et al. (2016). In principle, the higher the ED value is, the closer the object site. The raw ED was calculated at each SNP location using the equation:

A_aa, C_aa, T_aa, and G_aa represent the depth of bases A, C, T, and G on a site in the black grain bulk, respectively. A_ab, C_ab, T_ab, and G_ab represent the depth of bases A, C, T, and G on a site in the yellow grain bulk, respectively.

In order to increase the effect of large ED measurements and decrease the effects of background noise, the allele frequency of raw ED raised to the fifth power. Then the fitting result of ED⁵ calculated using local linear regression of the EDs with a span automatically chosen by minimizing the corrected Aikaike Information Criterion (AICc) (Hill et al., 2013), was used to associate analysis.

To minimize the genetic interval for fine-mapping and to verify the accuracy of SLAF-seq, We chose at about 1 Mb one to three potential SNPs located in the candidate region (Chr 1H 427,749,941 to 460,155,270 bp, IBSC 2012) to design their flanking primers using Oligo Primer Analysis Software v.7 which ranged from 100 to 300 bp in length. PCR reaction conditions were as follows: denaturation at 94°C for 5 min, 35 amplification cycles of 94°C for 30 s, annealing at 55°C for 30 s, and extension at 72°C for 30 s, with a final extension at 72°C for 5 min. PCR products were separated on 6% polyacrylamide gel (acrylamide/bisacrylamide ratio of 37.5:1) in 0.5 × Tris-Borate-EDTA (TBE) buffer and ran at room temperature for 2–4 h, stained with silver nitrate, and observed on white illumination. Size differences in polymorphisms were identified between Hatiexi No.1 and Zhe5819. PCR products showed no size polymorphisms on the polyacrylamide gel were sequenced in one direction using the specific PCR primers distal to the potential SNP position by biosune (Shanghai) Biotechnology Co., Ltd. The Megalign program (DNAStars) was used for sequence alignment and to confirm SNP sites.

The confirmed SNP markers were genotyped in 172 homozygous yellow individuals from F2 of Hatiexi No.1 × Zhe5819 following SNP marker detection with direct DNA sequencing or KASPar platform. Kompetitive Allele-Specific PCR (KASP) is a SNP genotyping system from LGC Genomics (United Kingdom) that tags different fluorescent dye to each SNP allele during the PCR reaction. Twenty two SNP markers were detected employing KASPar platform in segregating population by Beijing Vegetable Research Center (China). The KASP genotyping procedures were followed according Wen et al. (2015). The size differences markers were identified by polyacrylamide gel electrophoresis (PAGE). Other polymorphic markers were analyzed by Sanger DNA sequencing.

Linkage analysis of the molecular markers and black grain trait was performed using MAPMAKER version 3.0 software (Lander et al., 1987). Map distances were estimated using the Kosambi equation (Kosambi, 1944). For fine mapping, closer markers linked to the candidate gene were further developed and tested for their polymorphisms between the parents using Sanger DNA sequencing. Polymorphic markers were used for analysis of yellow grain plants from F2 generation. The alleles with the same genotype as that of black grain landrace Hatiexi No.1 was marked as ‘1,’ and alleles with the same genotype as that of yellow grain variety Zhe5819 was labeled as ‘0.’ For F2 plants with yellow grain, there are three possible genotypes for these markers, namely non-recombinants with ‘0/0,’ single recombinants with ‘0/1,’ and double recombinants with ‘1/1.’

After SLAF library construction and high-throughput sequencing, a total of 180,828,494 valid single-end reads were obtained, with each read length of ∼100 bp (Table 1). The GC content was 43.10% and the Q30 ratio was 92.92%. The SLAF numbers were 160,977 for Hatiexi No.1 and 181,313 for Zhe5819. The average sequence depths of SLAFs were ∼16.44- and ∼27.12-fold in black parent (Hatiexi No.1) and yellow parent (Zhe5819), respectively; and ∼45.41- and ∼41.27-fold in the black pool and yellow pool, respectively, (Table 1). SLAF tags were mapped on barley assembly (IBSC 2012) and 233,701 SLAFs markers distributed throughout the genomes. The SLAF numbers and chromosome positions are shown in Table 2.

From the 233,701 SLAF tags, 215,721 SNPs were obtained after aligning the sequence data to the barley reference. At the stage of SNP calling, SNPs with multiple allele loci and a depth less than 5× were filtered out. Polymorphic SNPs refer to SNPs that show polymorphic not only between the parents but also between the two bulked DNA samples. Finally, 77,542 polymorphic SNPs were ultimately selected for further analysis and the statistics of marker numbers on each chromosome according to the positioning result are shown in Table 2.

Both SNP_index and Euclidean distance association analysis were used to identify the candidate regions for barley black lemma and pericarp trait. For the SNP_index method, SNP_index was calculated for each identified SNP according to Abe et al. (2012) and Takagi et al. (2013). An average SNP_index of SNPs was calculated with 200 SNP_indexes located in a given genomic interval. SNP_index graphs were generated for the yellow (Figure 1A) and black (Figure 1B) pools by plotting the average SNP_index against the position of each sliding window in the barley genome assembly (IBSC 2012). After combining the SNP_index information into the yellow and black pools, the Δ(SNP_index) was calculated and plotted against the genome positions (Figure 1C). Peak regions above the threshold value were defined as those where Loess fitted values were greater than standard deviations above the genome-wide median in the Δ(SNP_index) plot. One candidate region associated with barley black grain spanned 49.28 Mb on chromosome 1H (from 414,847,463 to 464,122,721 bp, barley genome assembly, IBSC 2012), was identified with Δ(SNP_index) value above the threshold value of 0.26 (Figure 1C).

Euclidean distance (ED) was calculated for each SNP according Hill et al. (2013). To increase the effect of large ED measurements and decrease the effects of low ED measurements/noise, the 5th power of ED was calculated as the correlation value. The association threshold was 0.15 and one region on chromosome 1H was significantly correlated with the black lemma and pericarp trait. The result of the Euclidean distance association analysis was shown in Figure 1D. According to barley physical map (IBSC 2012), the candidate region was physically located on chromosome 1H between 427,749,941 and 460,155,270 bp, with a size of 32.41 Mb.

Combining the results of SNP_index and Euclidean distance association analysis suggested that the overlapped region (427,749,941–460,155,270 bp on chromosome 1H, IBSC 2012) was the candidate region of the barley black lemma and pericarp gene.

A total of 524 potential polymorphic SNPs were obtained in the 32.41 Mb candidate regions (Supplementary Table S1). To evaluate the accuracy of SLAF genotyping data, one to three SNPs per Mb were selected across the entire candidate region. Fifty four pairs of primers were designed due to their potential polymorphisms and physical position on barley genome assembly (IBSC 2012). Markers polymorphisms between Hatiexi No.1 and Zhe5819 were verified by electrophoresis and independent traditional Sanger sequencing. Thirteen of fifty four primer pairs showed no PCR products in one parent or both parents and were removed from analysis. HZSNP34 makers showed InDel polymorphism on polyacrylamide gels. The rest of the PCR products of the two parents showed no size polymorphisms on the polyacrylamide gels were sequenced directly. Sequences alignment between the parents identified twenty nine polymorphic markers (Supplementary Table S2). Among the twenty nine polymorphic markers, 24 markers showed SNP and five of them showed multi-nucleotide polymorphisms (Supplementary Table S2).

KASPar platform was used to conduct SNP genotyping in the F2 population consisting of 172 homozygous yellow grain individuals. Twenty two KASPar type SNP markers, including 19 SNP markers and 3 multi-nucleotide polymorphism markers (HZSNP15, HZSNP28, and HZSNP36), were designed (Supplementary Table S3). For three multi-nucleotide polymorphism markers, KASPar assays just screened one SNP, which is more than 50 bp away from the other variant sites. Except HZSNP28, all the KASPar type SNP markers genotyped the population successfully. InDel marker HZSNP34 was distinguished easily on 6% polyacrylamide gel in the population.

Linkage analysis showed that all markers were assigned to the target regions and the gene controlling black lemma and pericarp was delimited by markers HZSNP35 (1.9 cM) and HZSNP39 (6.2 cM) (Figure 2A). Moreover, the gene was co-segregated with HZSNP34 and HZSNP36. These results suggested that the markers mined from SLAF-seq-BSA data are reliable. According to the barley genome assembly (IBSC 2012), the markers order in the genetic map was not consistent with its physical map (Figure 2A and Supplementary Table S2). Thus, all marker sequences were blast against the current barley assembly released by the International Barley Sequencing Consortium in 2017 (IBSC 2017³). Blast alignment analysis showed that the genetic map was incompliance with the current physical map (IBSC 2017) and the physical distance between markers HZSNP35 (536444825–536445008 bp) and HZSNP39 (542121828–542122039 bp) to be approximately 5.68 Mb on IBSC 2017 assembly (Supplementary Table S2). The chromosomal location of this locus corresponded with black lemma and pericarp1 (Blp1) described by Costa et al. (2001) and Bungartz et al. (2016). Hence, we also named the gene as Blp1 following previously.

Markers developed by SLAF-seq in the 5.68 Mb (HZSNP35–HZSNP39) intervals were further screened to obtain polymorphic markers between the parents with direct DNA sequencing. Four polymorphic markers, including three co-dominant markers (HZSNP59, HZSNP61 and HZSNP63) and one dominant marker (HZSNP62) were identified (Supplementary Table S2). Then the co-dominant markers and HZSNP32 located in the reduced target region, were used to analyze the genotypes of yellow pericarp F2 plants. Among the 172 homozygous yellow F2 plants of Hatiexi No.1 × Zhe5819, six plants (Y225, Y314, Y333, Y372, Y401, and Y406) were recombinants on the HZSNP63 locus and nine plants (Y316, Y415, Y328, Y330, Y331, Y332, Y379, Y422, and Y444) on the HZSNP61 locus (Table 3). Two plants (Y316 and Y415) appeared to be recombinants on locus HZSNP59 in the downstream (Table 3). Because of the limited markers, no recombinant loci were found to be closer than HZSNP63. Eventually, the Blp1 gene was delimited within a 1.66 Mb (IBSC 2017 assembly, Chr 1H: 536,999,583-538,661,822) by the upstream marker HZSNP63 and the downstream marker HZSNP59 (Figure 2B and Table 3).