NGS sequencing and BSA-seq analysis

Total genomic DNA was extracted from bulked pools, and at least 3 μg genomic DNA was used to construct paired-end libraries with an insert size of 500 bp using the paired-end DNA sample prep kit (Illumina, San Diego, CA, USA). These libraries were sequenced using the HiSeq X10 (Illumina) NGS platform at Genedenovo (Guangzhou, China). Quality trimming is an essential step in generating high-confidence variant calling. Raw reads were processed to obtain high-quality clean reads according to three stringent filtering standards: 1) removing reads with ≥10% unidentified nucleotides; 2) removing reads with > 50% bases having phred quality scores ≤ 20; and 3) reads aligned to the barcode adapter.

To identify SNPs and indels, filtered reads were aligned to the Nipponbare reference genome sequence [53] using the Burrows–Wheeler Aligner (v.0.7.16a-r1181) with parameter ‘mem -M’; -M is an option used to mark shorter split-alignment hits as secondary alignments [54]. Variant calling was performed using the GATK UnifiedGenotyper (v.3.5; https://gatk.broadinstitute.org/hc/en-us/community/posts/360073637011-UnifiedGenotyper-in-GATK4). SNPs and indels were filtered using the GATK VariantFiltration function with proper standards (-Window 4, -filter "QD < 4.0 || FS > 60.0 || MQ < 40.0 ", -G_filter "GQ < 20"). All mutations were annotated for genes and functions, as well as genomic regions, using ANNOVAR [55]. Association analysis was performed using the SNP-Index [23], ∆(SNP-Index) [56], calculation of the G statistic [57, 58], ED [59], and two-tailed Fisher’s exact test [60] based on the SNPs. Average values were calculated using a 2000-kb sliding window with a step size of 20 kb and discarding of windows with less than 10 SNP/Indel. If the number of SNPs was insufficient, the results of this window were merged into the next window. The overlapping interval of the four methods was considered as the final QTL interval of CT.