Forty-six individuals, derived from the tails of a bi-parental recombinant inbred line population (RIL_F6:7), previously studied for five high-correlated morpho-physiological traits (i.e., plant height, juvenile growth habit, total tiller number, fertile tiller number and heading date) (Vitale et al., 2021), were analyzed. The population of 176 RILs was developed from a cross between Lankaodali and Rebelde bread wheat cultivars, using the single seed descent method by advancing random F2 plants to the F6:F7 generation (Vitale et al., 2021). The two parental accessions were characterized by contrasting quantitative and qualitative agronomic traits (Botticella et al., 2018; Vitale et al., 2021).

Lankaodali is an early bread wheat cultivar of Chinese origin with very large kernels, low tillering ability, and poor qualitative attributes. Rebelde is an Italian cultivar with late flowering, small kernels, high tillering ability and excellent grain quality traits. During the 2020–2021 growing season, the leaves of 14-day-old seedlings of each of 48 individuals were collected and ground using liquid nitrogen. DNA was extracted using the Quick-DNA Plant/Seed Miniprep Kit (Zymo Research, United States) according to the manufacturer’s instructions. DNA quality and quantity were estimated using the NanoDrop ND-1000 spectrophotometer (Thermo Scientific, Wilmington, DE, United States) and the Qubit fluorometer (Invitrogen, Carlsbad, CA, United States), respectively.

In-solution-based hybridization was applied to capture the target loci. Baits (304,327) were designed to specifically capture over 250 Mb of coding DNA sequences (CDS) (myBaits®, Arbor Biosciences, Ann Arbour, MI, United States; http://www.arborbiosci.com). Bait design was based on the Chinese spring wheat genome v1.0 (IWGSC et al., 2018) and was carried out following the manufacturer’s protocol version 3.01 (https://arborbiosci.com/wp-content/uploads/2020/01/myBaitsExpert_WheatExome_Product_Sheet.pdf; accessed 14/07/2022).

As shown in Supplementary Table S1, the baits were evenly distributed along the genome, with the lowest number (10,345 baits) on chromosome 4D and the highest (17,587 baits) on chromosome 2B.

The extracted DNA (250 ng −1 μg) was subjected to random mechanical shearing to obtain fragments with an average size of 400 bp. The fragments underwent an A-tailing reaction at 3’ of the blunt-end, where barcoded adapters were then ligated. Libraries were paired-end sequenced on an Illumina NovaSeq 6,000 platform. After Illumina sequencing, an average of 153 million raw reads per individual were obtained, ranging from 96 M (sample 21) to 230 M reads (sample 31) (Supplementary Table S2). Overall sequencing quality was assessed by FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and reads with a base quality score below 20 (Q < 20) were removed using Trimmomatic (Bolger et al., 2014) (SLIDINGWINDOW = 6). High quality reads were then aligned to the Chinese spring wheat genome v1.0 (IWGSC et al., 2018) with BWA-MEM (Li et al., 2009) (MINIMUMSEEDLENGTH = 19; BANDWIDTH = 100). SAMtools (Heng et al., 2009) were used to convert SAM files to BAM; the latter were processed by the MarkDuplicates utility of Picard version 1.109, (http://picard.sourceforge.net), to remove duplicate reads. Reads with mapping quality scores below 30 (Q < 30) were also filtered out. QualiMap (Garcia Alcalde et al., 2009) was used to evaluate the effectiveness of the mapping process: the percentage of mapped reads was greater than 99% for each individual with an average error rate of 0.74% (Supplementary Table S2). The reads-to-genome mapping produced over 55% of on-target reads for all individuals except sample 26, which showed the lowest percentage (∼47%) (Supplementary Figure S1). Hybridization capture methods are known to be prone to off-target enrichment and capture (Kaur et al., 2017) and in polyploid species target specificity and efficiency is affected by the presence of homoeologous sequences (King et al., 2015). This reflects the percentage of on-target reads we have achieved which is consistent with what obtained by King et al. (2015).

The CoverageBed utility of the BEDtools package (http://bedtools.readthedocs.org/) was used to derive the depth of coverage of the target regions for each individual. On average, more than 50% of the bases in the bait regions were covered at a depth greater than 20x (Figure 1A). Coverage per individual averaged 13.7x, with most individuals showing a mean coverage between 4x and 35x (Supplementary Table S3). We also estimated the depth of coverage per-base at the gene level (Figure 1B). As an example, Figure 1B shows per-base depth of coverage for the TraesCS2B01G175300 gene in all individuals. As expected, the depth of coverage was fairly uniform across exons and very shallow for introns.

SNPs were identified using the Genome Analysis Toolkit (GATK) version 4.0 (Van der Auwera et al., 2013), following the best practices recommended by the GATK documentation. The VariantRecalibrator utility was used to distinguish true genetic variants from sequencing or data processing artifacts (false positives) and remove the latter from the raw VCF (Variant Call Format) file. A total of 15,046,465 SNPs with a sequencing and alignment quality score ≥30 (Phred-scaled) and a coverage ≥10 were identified. The VCF-stat utility in VCFtools (King et al., 2015) was used to retrieve the statistics by sample. The number of variants ranged from 509,707 SNPs (sample 21) to 1,939,789 (sample 37), with a mean of 1,062,065 SNPs. To assess the level of residual heterozygosity for each individual, we calculated the percentage of heterozygosity at each locus by dividing the number of heterozygous genotypes for a given locus by the total number of individuals. This number was averaged across all the chromosomes and returned residual heterozygosity equal to 5.85% (Supplementary Figure S2) ( Danecek et al., 2011). As expected, transitions (Ts) were the most abundant (∼70%), while transversions (Tv) accounted for ∼30%. Their ratio was 2.11. SnpEff (Cingolani et al., 2012; Velmurugan et al., 2018) was used to classify variants based on genomic location and biological effect (Supplementary Figures S3, S4). As observed by Suren et al. (2016), we found that some of the variants also came from non-target regions, although bait design was exclusively based on annotated exons. Indeed, most SNPs were detected in the intergenic (47%) and upstream gene regions (17%), while approximately 26% fell within exons (Supplementary Figure S3). A stacked bar chart showing SNPs grouped by their impact (low, moderate, and high) is reported in Supplementary Figure S4. SNPs with “moderate” impact on protein functioning were the most abundant (52%). Variants marked as “low” impact accounted, on average, for 46% of all variants (Supplementary Figure S4). Finally, the variants estimated to have a deleterious impact on gene functioning (i.e., “high” impact) accounted for 2% on average. The latter are to be considered the most interesting variants, since they might determine phenotypes of interest.

The conversion of SNP markers into KASP markers is particularly challenging in polyploid crops, due to the presence of homoeologous sequences (Makhoul et al., 2020). To validate the quality of polymorphism discovery and variant annotation generated by exome capture, we designed and tested eight KASP assays that distinguished between alleles for the 48 individuals (Figure 2). SNP markers were selected following different criteria: presence of polymorphism between parental lines, sequencing depth, biological effect, randomness. The KASP primers design was carried out using the commercial KASP assay design service (KASP-by-Design) developed by LGC Biosearch Technologies. For each SNP, 100-bp flanking regions were obtained from the reference genome and converted into KASP to detect the specific parental allele. The LGC Genomics (Hoddeson, United Kingdom) designed two allele-specific forward primers carrying the standard FAM- or VIC-compatible tail (FAM: 5′-GAA​GGT​GAC​CAA​GTT​CAT​GCT-3′; VIC: 5′-GAA​GGT​CGG​AGT​CAA​CGG​ATT-3′) with the targeted SNP at the 3′ end (Supplementary Table S4). Almost all primer pairs had perfect match with the target sequences and mismatches at 3’ tail with the other wheat sub-genomes. DNA samples were arrayed into a 96-well PCR plate, each containing ∼5 µl reaction mix (45 ng of dry DNA, 2.5 µl of 1 × KASP master mixture, and 0.1 µl of primer mix). Primer mix included a final concentration of 30 µM of the common primer and 10 µM of each tail primer. PCR experiments were performed using the ABI ViiA7 instrument (Applied Biosystems, Foster City, CA, United States) as follows: initial denaturation at 94°C for 15 min, followed by ten touchdown cycles (94°C for 20 s; touchdown at an initial temperature of 61°C and decreasing by −0, 6°C per cycle for 60 s) reaching a final annealing temperature of 55°C, followed by 26 additional annealing cycles (94°C for 20 s; 55°C for 60 s). Fluorescence readings were performed with a temperature below 37 °C and allelic discrimination plots were drawn using the SNP viewer software (https://www.biosearchtech.com/support/tools/genotyping-software/snpviewer). Individuals with contrasting alleles at each SNP locus (0/0 corresponds to the reference allele, 1/1 to the alternative allele, 0/1 to a heterozygous locus, and “N/A” to a missing data point) were genotyped using the KASP assays (Supplementary Table S5). Results showed that roughly 80% of alleles for all loci were scored identically between WES data and KASP genotyping, including those with low coverage (Supplementary Table S5).

In addition to sequence data and sequence variations, we will make available the plant material described in this paper under a standard material transfer agreement This will allow the scientific community to identify marker-trait associations linked to key agronomic traits (including those mentioned above) assessed in other environments (year, site, management), and, thus, detect stable QTL based on genotype-environment (GxE) interactions.

Indeed, SNP discovery supports the use of bulked segregant analysis as a powerful tool to accelerate gene identification and QTL mapping cost-effectively.

In addition, the data described here could be used to dissect naturally occurring allelic variation at candidate genes controlling key agronomic traits, identify useful alleles (i.e., loci affecting the traits of interest) and understanding their function.