Benchmarking of low coverage sequencing workflows for precision genotyping in eggplant

RESEARCH Open Access © The Author(s) 2025. Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit h t t p : / / c r e a t i v e c o m m o n s . o r g / l i c e n s e s / b y - n c - n d / 4 . 0 /. Baraja-Fonseca et al. BMC Plant Biology (2025) 25:1125 https://doi.org/10.1186/s12870-025-07242-x BMC Plant Biology *Correspondence: Virginia Baraja-Fonseca [email protected] Pietro Gramazio [email protected] 1Instituto de Conservación y Mejora de la Agrodiversidad Valenciana, Universitat Politècnica de València, Camino de Vera 14, Valencia 46022, Spain 2Instituto de Biología Molecular y Celular de Plantas, Consejo Superior de Investigaciones Científicas - Universitat Politècnica de València, Camino de Vera 14, Valencia 46022, Spain Abstract Background Low-coverage whole-genome sequencing (lcWGS) presents a cost-effective solution for genotyping, particularly in applications requiring high marker density and reduced costs. In this study, we evaluated lcWGS for eggplant genotyping using eight founder accessions from the first eggplant MAGIC population (MEGGIC). We tested various sequencing coverages and minimum depth of coverage thresholds with two SNP callers, Freebayes and GATK. Reference SNP panels were used to estimate the percentage of common biallelic SNPs (i.e., true positives) relative to the low coverage datasets (accuracy) and the SNP panels themselves (sensitivity). Furthermore, the percentage of true positives with the same genotype across both datasets was calculated to assess genotypic concordance. Results Sequencing coverages as low as 1X and 2X achieved high accuracy but lacked sufficient sensitivity and genotypic concordance. However, 3X sequencing reached approximately 10% less sensitivity than 5X while maintaining genotypic concordance above 90% at any depth of coverage threshold. Freebayes outperformed GATK in terms of sensitivity and genotypic concordance. Therefore, we used this software to conduct a pilot test with some MEGGIC lines from the fifth generation of selfing, comparing their datasets with a gold standard. Sequencing coverages as low as 1X identified a substantial number of true positives, with 3X significantly increasing the yield, particularly at moderate depth of coverage thresholds. Additionally, at least 30% of the true positives were consistently genotyped in all lines when using coverages greater than 2X, regardless of the depth of coverage threshold applied. Conclusions This study highlights the importance of using a gold standard to reduce false positives and demonstrates that lcWGS, with proper filtering, is a valuable alternative to high-coverage sequencing for eggplant genotyping, with potential applications to other crops. Keywords Eggplant (Solanum melongena), Genotyping, Low-coverage whole-genome sequencing (lcWGS), Bioinformatic pipeline, Benchmarking analysis, Gold standard (GS) Benchmarking of low coverage sequencing workflows for precision genotyping in eggplant VirginiaBaraja-Fonseca1*, AndreaArrones1, SantiagoVilanova1, MariolaPlazas1, JaimeProhens1, AurelianoBombarely2 and PietroGramazio1* Page 2 of 16Baraja-Fonseca et al. BMC Plant Biology (2025) 25:1125 Background Plant genomic characterization is a critical step in modern breeding programs, and essential for studying diversity, understanding domestication and recombination events, and identifying candidate regions linked to important agronomic traits [1]. Current methods for high-throughput genotyping in plants primarily involve two approaches: reduced representation sequencing (RRS) and whole-genome sequencing (WGS) [2]. However, low-coverage whole-genome sequencing (lcWGS) is gaining popularity as a genotyping strategy that combines the broad genomic coverage of WGS with the cost-efficiency of RRS [3, 4]. lcWGS allows for the identification of a dense marker set with a comprehensive representation of the entire genome, using very low sequencing coverage (< 10X) [5, 6]. This strategy has opened new possibilities for genotyping studies in various crops, including chickpea [7], rice [8], canola [5], soybean [6], tomato [9], radish [10], wheat [11] and potato [12], among others. The primary strength of lcWGS lies in its cost-effectiveness, with costs scaling down in tandem with reduced sequencing coverage [4, 13]. Another valuable feature of this strategy is the concurrent reduction in data volume as sequencing coverage diminishes, expediting and streamlining bioinformatic analysis [5, 14]. Nevertheless, low sequencing coverages can lead to erroneous conclusions due to the limited information provided by the low number of reads. Some potential challenges encompass: (1) genotype misclassification, (2) loss of genuine polymorphism, and (3) sequencing errors being erroneously classified as genetic variants [15]. To address these drawbacks, rigorous SNP filtering steps are critical [13, 16]. Furthermore, employing additional procedures is advised for the elimination of false positive calls, such as using more than one SNP calling software [17, 18] and validating polymorphisms against a set of truly-assumed genetic variants (gold standard; GS), supported by a higher number of reads or validated in various independent studies [19–21]. Eggplant (Solanum melongena L.) is an economically significant crop ranking as the third most important solanaceous crop after potato and tomato in global production and the fifth among all vegetable crops [22]. Despite its economic importance, available genetic and genomic resources of eggplant have traditionally lagged behind those of other important vegetable crops [23]. However, noteworthy progress includes the development of the first and only multiparent advanced generation inter-cross (MAGIC) population in eggplant, known as MEGGIC [24]. To fully leverage MAGIC populations as valuable next-generation genomic resources, genotypic characterization is essential. The 5k Single Primer Enrichment Technology (SPET) genotyping platform [25], developed from the resequencing at 20X of its eight founders [26], was used to genotype 420 individuals from the third generation of selfing (S3MEGGIC), resulting in 7724 high-confidence SNPs [24]. Even though the SPET genotyping allowed the dissection of key genes for eggplant genetics and breeding [24, 27], the genetic characterization of the segregating individuals through variant identification and haplotype resolution was not fully comprehensive. Two principal limitations of this genotyping approach by amplicon sequencing are the limited number of variants that can be interrogated and their distribution, which are preferentially selected in the gene-rich chromosome arms to assess gene allelic diversity [25]. These issues can be addressed by lcWGS, as demonstrated in recent specific eggplant studies that created high-density recombination bin-based genetic maps and improved QTL mapping resolution [28, 29]. However, there remains a limited understanding of the impact of several parameters during data processing on its accuracy and sensitivity. Thus, the principal aim of this study is to establish an optimized workflow for analysing low-coverage genomic data in eggplant. We used the MEGGIC founders to benchmark different combinations of sequencing coverages, minimum depth of coverage (DP) thresholds and SNP callers. A proof-of-concept validation of these findings was conducted on lines from the fifth generation of selfing of the eggplant MAGIC population (S5MEGGIC), which will facilitate the optimization of genomic diversity analyses in eggplant collections and populations. Furthermore, this study provides guidelines for selecting appropriate parameters in eggplant genomics analysis and presents a protocol that can be broadly applied across various crops and research objectives. Methods Plant materials, library preparation and resequencing The plant materials used for the low-coverage sequencing benchmarking were the eight founders of the MEGGIC population, consisting of seven Solanum melongena and one wild relative S. incanum accessions (Additional file 1) [24]. All the accessions are maintained at the Universitat Politècnica de València (UPV) germplasm bank. On the other hand, to increase the genomic characterization precision of the MEGGIC population and to validate the lcWGS benchmarking results of this study, we used four random recombinant lines from the S5 generation (labelled for this study as S5-1, S5-2, S5-3, S5-4). The S5 lines were obtained following a funnel scheme as described by Mangino et al. [24] (Additional file 1). Seeds from the twelve samples (the eight founders and the four S5MEGGIC lines) were germinated in Petri dishes, following the protocol developed by Ranil et al. [30]. We then transferred them to seedling trays in Page 3 of 16Baraja-Fonseca et al. BMC Plant Biology (2025) 25:1125 a climatic chamber under a photoperiod and temperature regime of 16h light (25°C, 100–112 µmol m− 2s− 1) and 8h dark (18°C). Total genomic DNA was extracted from approximately 100 mg of young leaves following the SILEX protocol described by Vilanova et al. [31]. We assessed DNA integrity and quality of the extracted DNA through agarose electrophoresis and NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, Delaware, USA). DNA concentration was determined with a Qubit® 2.0 fluorometer (Thermo Fisher Scientific, Waltham, MA, USA). High-quality DNA samples (260/280 and 260/230 ratios > 1.8) were then shipped to the Beijing Genomics Institute (BGI Genomics, Hong Kong, China) for the construction of 150bp paired-end libraries and subsequent sequencing using the DNBseq platform. The lcWGS was pursued to generate around 6.3 Gb high-quality sequence reads per sample, approximating 5X genome coverage. Raw reads underwent filtration using SOAPnuke software [32] to remove adapters and low-quality reads (-n 0.001 -l 10 -q 0.4 --adaMR 0.25 --ada_trim). After receiving trimmed reads, we performed quality control using FastQC (version 0.11.9) [33] to assess the effectiveness of the quality filtering (Fig.1.A and Fig.1.B). Downsampling and mapping The original fastq files obtained from the 5X sequencing were utilized to simulate different lower-depth samples using seqtk tool (version 1.3-r106, h t t p s : / / g i t h u b . c o m / l h 3 / s e q t k ) . Each sample was computationally subsetted based on the number of reads. Average depths of 1X, 2X, 3X and 4X were produced (Fig.1A and B). The same random seed (-s) was employed to preserve read pairing. For each simulated dataset, five replicates (R) were randomly generated (Fig.1A and B). We mapped clean reads from the five gradient sequencing coverages for each sample against the v3.0 “67/3” high-quality eggplant reference genome [34] using BWA with its Minimal Exact Match algorithm (BWA-MEM) (version v.0.7.17–r1188; Fig.1A and B) [35]. The resulting alignment data were subsequently transformed into the BAM format using the SAMtools package (version 1.13) [36]. Mapping statistics, including the number of mapped and unmapped reads, and the average depth of coverage, were recorded in the output file generated by the QualiMap application (version 2.2.1) [37]. Moreover, we used the ‘coverage’ function within the SAMtools package (version 1.13) [38] to calculate genome coverage. In-depth analysis of the spatial distribution of mapping depth across the genome, with a window size of 10 Kbp, was performed with the bamCoverage tool (version 3.5.1) [39]. Output data was visualized and graphically represented using the ‘plot’ function (version 3.6.2) in R (version 4.3.2) [40]. Finally, we marked PCR duplicates using the MarkDuplicates tool from Picard software (version 1.119; h t t p s : / / b r o a d i n s t i t u t e . g i t h u b . i o / p i c a r d /; Fig. 1A and B). Polymorphism detection and data filtering Variant calling in low-coverage founders’ data was carried out at the sample level (i.e., each sample independently for each sequencing coverage) using two different Bayesian-based software: Freebayes (version 1.3.6) [41] and GATK HaplotypeCaller (version 4.3.0.0) (Fig. 1A) [42]. BAM files from the mapping step were provided as input to both tools individually, using the -b and -I options, respectively. We ran both callers with default configuration settings, except for the minimum quality requirements for mapping and base, for which we set the threshold at 20. Biallelic SNPs, excluding monomorphic ones, were kept using BCFtools (version 1.13; h t t p s : / / s a m t o o l s . g i t h u b . i o / b c f t o o l s / b c f t o o l s . h t m l). To assess the impact of the minimum depth of coverage thresholds on polymorphism detection, we tested a range of thresholds from DP 1 to DP 10. This range was selected to capture sufficient genomic variation while balancing the need for accuracy in identifying true positives and minimizing false positives. This process generated 400 list-based SNP sets per low-coverage level-SNP caller combination. Specifically, this includes eight samples, five replicates, and 10 different minimum depth of coverage thresholds. For the 5X coverage level, 80 SNP sets were generated per SNP caller, as no replicates were included for this coverage (Fig.1A). Variants of the four S5MEGGIC lines were identified at the population level (i.e., the four samples together) using Freebayes (Fig. 1B). BAM files from the mapping step were provided as a list of inputs, using the -L option. This process created a single output file containing information for all the samples. While we maintained the default configuration settings, we specifically set the thresholds for mapping and base quality to a minimum of 20. In the same way, as did with founders’ data, biallelic SNPs were filtered by the minimum depth of coverage ranging from DP 1 to DP 10. Finally, only polymorphic variants among the four lines were retained (Fig.1B). MEGGIC founders derived reference SNP panels and gold standard To benchmark each MEGGIC founder’s combinations of lc-DP (low-coverage dataset and minimum depth of coverage threshold), two reference SNP panels were established for each MEGGIC founder from its 20X resequencing dataset (SRA BioProject PRJNA392603) (Fig. 1C) [26]. The cleaned reads underwent a quality control assessment using FastQC (version 0.11.9) [33]. Then, we conducted the mapping and polymorphism detection at the sample level using Freebayes (Freebayes Page 4 of 16Baraja-Fonseca et al. BMC Plant Biology (2025) 25:1125 SNP panel) and GATK (GATK SNP panel) as previously detailed for the low-coverage founders’ data. Biallelic SNPs supported by at least 20 reads were retained to ensure the most informative markers, while monomorphic ones were excluded (Fig.1C). In addition, we established a unique gold standard to select common biallelic SNPs identified in the S5MEGGIC lines. This process validated the benchmark results and determined the optimal combination of parameters for the genomic characterization of the S5MEGGIC Fig. 1 Bioinformatics pipeline for low-coverage genomic data analysis. Squares represent files and circles indicate procedures during the bioinformatic analysis. A Low-coverage whole genome sequencing benchmark workflow. MEGGIC founders’ 5X cleaned data were down-sampled to 1X, 2X, 3X and 4X. Five replicates were generated for each down-sampled level (R1-R5). After the mapping step, SNP calling was performed for each low-coverage dataset and founder using Freebayes (F) and GATK (G), followed by variant filtration based on the minimum depth of coverage (DP 1 to DP 10). The SNPs datasets were validated using the corresponding reference SNP panel. B Proof-of-concept validation of the benchmark study. S5MEGGIC lines’ 5X cleaned data were downsampled to 1X-4X, and five replicates were generated for each level (R1-R5). After the mapping step, Freebayes was used to perform population-level SNP calling on the combined BAM files from all lines for each sequencing coverage and replicate. Variant filtration was based on the minimum depth of coverage (DP 1 to DP 10). The SNP datasets were validated using the gold standard (GS). Output datasets were labelled according to the sequencing coverage and the applied filter (e.g., 1X DP 1 indicates a sequencing coverage of 1X and a minimum depth of coverage threshold set at 1). C Reference SNP panels and GS preparation from MEGGIC founders’ 20X data. One reference SNP panel per genotype-SNP caller (Freebayes, F; GATK, G) combination was obtained. Additionally, a unique GS was obtained by performing a population-level SNP calling using Freebayes followed by variant filtration Page 5 of 16Baraja-Fonseca et al. BMC Plant Biology (2025) 25:1125 population (Fig.1C). The GS was derived by performing SNP calling at the population level, as previously described for the low-coverage S5MEGGIC lines’ data, using the 20X resequencing datasets from all eight founders with Freebayes (-q 20 -m 20 –limit-coverage 800). The resulting VCF file included both common and unique variants across all founder genomes. Biallelic SNPs supported by at least 20 reads were retained, and monomorphic sites among samples were removed (Fig.1C). Genotyping accuracy, sensitivity and concordance evaluation To determine the optimal combination of tools and parameters in terms of variant discovery, we compared the low-coverage datasets (lc-datasets) from each founder to their respective Freebayes and GATK SNP panels, with the latter serving as the reference due to their superior read support (Fig. 1A). These comparisons were performed with the isec tool from the BCFtools package (version 1.13; h t t p s : / / s a m t o o l s . g i t h u b . i o / b c f t o o l s / b c f t o o l s . h t m l) with the default configuration, yielding only records with identical alleles between the reference SNP panel and the lc-dataset (bcftools isec -c none). The output files allowed the identification of polymorphic sites: (I) common between the reference SNP panel and the lcdataset (true positives; TP), (II) private to the lc-dataset (PlcD), and (III) private to the reference SNP panel (PP). The isec tool was also used to accurately determine the common variants (TP) between the S5MEGGIC lines and the GS, assuming the genotypes in the full 20X data were correct due to the high number of reads supporting each variant (Fig.1B). The first metric evaluated was accuracy, defined as the capability to filter out correctly potential false polymorphisms or identification errors. It was calculated as the ratio of true positives to the total number of polymorphisms identified in the sample (1). Accuracy = TP TP +PlcD (1) On the other hand, sensitivity was assessed as the capability to identify correctly genuine polymorphisms within the genome. It was calculated as the ratio of true positives to the total number of polymorphisms identified in the reference SNP panels (2). Sensitivity = TP TP +PP (2) Finally, genotypic concordance was calculated as the percentage of true positives in the lc-dataset that exhibited complete allele matches (homozygous-homozygous or heterozygous-heterozygous) with the genotype at the same site sequenced at 20X coverage. This measure reflected the ability to accurately assign genotypes even at low coverages. Statistical analysis Multifactorial analysis of variance (ANOVA) was used to assess the effects of SNP caller, coverage level, and minimum depth of coverage threshold on accuracy, sensitivity and genotypic concordance. Then, we conducted posthoc pairwise comparisons using Tukey’s Honest Significant Difference (HSD) test to identify specific differences between the levels of each factor. The significance level was set at p < 0.05. All analyses were carried out using Statgraphics Centurion 19 software (Statgraphics Technologies, Inc., The Plains, VA, USA). Results LcWGS and mapping The 5X lcWGS of the seven S. melongena and one S. incanum founders of the MEGGIC population, along with four recombinant S5MEGGIC lines, yielded an average of 19.70M reads per sample (Additional file 2). Of these, 96.38% were high-quality nucleotide bases, with a Q score greater than 20 (> Q20) (Additional file 2). To benchmark the impact of varying sequencing coverages, subsets ranging from 1X to 4X were simulated (Fig.1A and B). The original 5X and low-coverage datasets were aligned against the “67/3” eggplant reference genome using BWA-MEM software. The average percentage of mapped reads was 96.40%, with no significant differences observed between the relative data of the 5X and lowcoverage sets (Additional file 3). The BAM files derived from the original 5X sequences exhibited an average depth of coverage of 4.79, and the 1X, 2X, 3X, and 4X subsets of 0.96, 1.92, 2.87, and 3.83, respectively (Additional file 3). As the average sequencing coverage increased, the proportion of the reference genome covered also expanded, although the increments became progressively smaller at higher coverages (Fig.2A and Additional file 3). At 1X, 39.06% of the reference genome was covered, increasing to 55.20% at 2X and 63.05% at 3X. However, the gains diminished significantly at higher coverages, with only 66.54% and 68.52% covered at 4X and 5X, respectively (Fig. 2A). Surprisingly, a marginal increment of 4.58% was observed from 5X to 20X coverage (Fig.2A). The distribution of mapped reads across the sequencing coverages displayed a consistent pattern, characterized by regions of the genome lacking sequencing data (underserved regions) alongside areas subjected to excessive sequencing (over-sequenced regions) (Fig. 2B and Additional file 4). Page 6 of 16Baraja-Fonseca et al. BMC Plant Biology (2025) 25:1125 In silico evaluation of LcWGR under varying parameters To assess the effectiveness of lcWGS for polymorphism identification, we systematically tested several combinations of sequencing coverages (1X to 5X), depth thresholds (DP 1 to DP 10), and two widely-used SNP callers, Freebayes and GATK (Additional file 5). These comparisons aimed to determine which SNP caller and parameters combination provided the best balance between polymorphism yield and resource efficiency. Regarding SNP caller performance, five key trends emerged from our analysis: (I) Freebayes identified more polymorphisms than GATK, particularly at higher sequencing coverages and lower DP thresholds (Additional file 6 and Additional file 7); (II) while Freebayes showed no significant differences between filtering at DP 1 and DP 2, GATK exhibited a slight reduction in biallelic SNPs, particularly at lower sequencing coverages (Fig.3A and Additional file 7); (III) the reduction in the number of SNPs identified from DP 1 to DP 10 was more Fig. 3 A Variation in the percentage of total biallelic SNPs when using different minimum depth of coverage thresholds (from DP 2 to DP 10) compared to DP 1 (set at 100%) for different sequencing coverages (1-5X) using Freebayes and GATK. Percentages represent the average across the eight founder accessions, with SD indicated by error bars (n = 40 for 1-4X and n = 8 for 5X). B Variation in the scaling factor, indicating the change in the number of total biallelic SNPs when increasing the sequencing coverage from one coverage to the next, using different minimum depth of coverage thresholds (from DP 1 to DP 10) with Freebayes and GATK. Values represent the average across the eight founder accessions, with SD indicated by error bars (n = 40) Fig. 2 A Percentage of reference genome coverage at different sequencing coverages (1X to 5X, and 20X) and the increments (yellow bars) from one coverage to the next. Genome coverage refers to the percentage of the genome covered by at least one read. The data shown are averages calculated from five replicates of each genotype for 1-4X sequencing coverage and from the original datasets for each genotype for 5X and 20X coverage. Standard deviation is indicated by error bars (n = 60 for 1-4X, n = 12 for 5X and 20X). B Example of the distribution of mapped read across the parental A chromosome 1 by sequencing coverage. The maximum cutoff is set at DP 20 for low coverages and DP 60 for 20X data. Peaks represent regions of high sequencing within 10 Kbp windows Page 7 of 16Baraja-Fonseca et al. BMC Plant Biology (2025) 25:1125 Fig. 4 (See legend on next page.) Page 8 of 16Baraja-Fonseca et al. BMC Plant Biology (2025) 25:1125 pronounced in GATK (88.92%) compared to Freebayes (81.63%), indicating a more aggressive filtering effect in GATK at higher DP thresholds (Fig.3A); (IV) the depth of coverage plateau was reached earlier with GATK than with Freebayes (e.g., for 1X, the difference between filtering at DP 6 and DP 10 was 5.60% for GATK and 10.91% for Freebayes) (Fig. 3A), which aligns with the lower number of reads used by GATK (Additional file 8); and (V) GATK exhibited a higher scaling factor between sequencing coverages, meaning it identified a larger proportion of biallelic SNPs at each coverage level relative to the previous one (e.g., 2X:1X). Specifically, GATK identified a larger proportion of new SNPs as coverage increased from 2X to 3 × (2.16) and 4X to 5 × (1.45), compared to Freebayes (1.97 and 1.34, respectively) (Fig.3B). However, Freebayes slightly surpassed GATK when moving from 1X to 2 × (3.31 vs. 3.29) (Fig.3B). Both sequencing coverage and depth of coverage threshold had a significant influence on the number of polymorphisms identified, with higher coverage and lower thresholds consistently yielding more biallelic SNPs (Additional file 7). However, the impact of the depth of coverage threshold was dependent on the sequencing coverage. Specifically, the differences could be attributed to how SNPs were distributed across different depth levels (Additional file 8). At lower coverage, the SNP reduction was more pronounced at lower DP thresholds than at higher ones (Fig.3A). For example, at 1X with Freebayes (median = DP 3, and 3rd quartile ≤ DP 4), the SNP reduction from DP 2 to DP 3 was 38.54%, whereas the difference between DP 9 and DP 10 was only 1.15% (Fig.3A). Conversely, at higher sequencing coverages, the impact of DP thresholds on SNP reduction was less significant. For instance, at 5X with Freebayes (median = DP 8, and 3rd quartile ≤ DP 11), filtering at DP 3 resulted in a 5.71% reduction compared to DP 2, which was closer to the 7.48% difference observed between filtering at DP 9 and DP 10 (Fig.3A). This trend was also observed in the scaling factor when transitioning from one sequencing coverage level to the next. When comparing the number of SNPs identified at 1X and 2 × (2X:1X), the scaling factor was around 2 at DP 1 and DP 2, indicating that at 2X, the number of SNPs identified was double that at 1X (Fig.3B). At higher depth of coverage thresholds, the scaling factor increased to 3 or 4 with both Freebayes and GATK (Fig.3B). In contrast, the scaling factor was more consistent across DP thresholds at higher sequencing coverages. For example, when moving from 4X to 5X, the scaling factor ranged from 1.16 at DP 1 to 1.59 at DP 10 using Freebayes, and from 1.18 to 1.80 with GATK (Fig.3B). Additionally, at lower sequencing coverages, a plateau was achieved for the scaling factor, but this was not observed at higher coverages (5X:4X), where the scaling factor increased steadily across depth of coverage thresholds. Prediction of the optimal combination of factors to perform LcWGS Comparisons between the reference SNP panels and each SNP dataset across lc-DP combinations allowed for the determination of true positives, missing data, and genotype assignment errors, as well as the estimation of accuracy and sensitivity. These panels were generated from the 20X resequencing dataset of the MEGGIC founders [26], assuming that the genotypes were more accurately characterized due to the higher coverage. The total cohorts of biallelic SNPs constituting the Freebayes SNP panels ranged from 3.68M to 7.59M (Additional file 9). On average, these were 1.58 times greater than those identified by GATK, which ranged from 1.86M to 4.93M (Additional file 9). Per sample-SNP caller combination, we evaluated 210 SNP datasets generated from 5X downsampling, observing a consistent increase in true positives with higher sequencing coverage and a reduction with stricter depth of coverage thresholds (Additional file 7). Freebayes consistently identified more true positives than GATK across all coverage levels, which aligned with results observed at 20X (Additional file 9). For instance, at DP 1, Freebayes outperformed GATK, identifying an additional 231.20 k TP at 1X, 574.71 k at 3X and 3.26M at 5X. Similarly, at the more stringent DP 10 threshold, Freebayes detected more TP than GATK across coverages, with differences ranging from 11.20 k at 1X to 425.03 k at 5X (Additional file 7). Accuracy trend differed between the two SNP callers (Additional file 10). At low depth of coverage thresholds (DP 1 to DP 4), Freebayes achieved higher average accuracy than GATK. While Freebayes exhibited values that ranged from 44.82% at 5X DP 1 to 49.87% at 1X DP 1, the average accuracy obtained with GATK at the same depth of coverage threshold showed less variability among sequencing coverages, ranged from 39.07% at 5X to 39.41% at 1X (Fig.4.A). However, at higher depth of (See figure on previous page.) Fig. 4 Accuracy, sensitivity, and genotypic concordance and discordance achieved for each combination of sequencing coverage (1-5X) and minimum depth of coverage threshold (from DP 1 to DP 10) using Freebayes and GATK. Percentages represent the average across the eight founder accessions, with SD indicated by error bars (n = 40 for 1-4X and n = 8 for 5X). A Accuracy refers to the ability to correctly filter out potential false positives or identification errors, calculated as the ratio of true positives to the total number of identified polymorphisms. B Sensitivity represents the ability to detect genuine polymorphisms within the genome, calculated as the ratio of true positives to the total number of polymorphisms identified in reference SNP panels. C Percentage of true positives with identical genotypes between the lc-datasets and the reference SNP panels. (D) Percentage of heterozygous true positives misclassified as homozygous in the lc-datasets. E Percentage of homozygous true positives misclassified as heterozygous in the lc-datasets Page 9 of 16Baraja-Fonseca et al. BMC Plant Biology (2025) 25:1125 coverage thresholds, GATK not only surpassed Freebayes in accuracy but also exhibited increased accuracy variability among sequencing coverages (Additional file 10). This trend was evident from DP 6 at 1X, continuing to increase up to DP 10 for 5X coverage (Fig.4). Notably, Freebayes achieved its highest values (> 59%) at specific thresholds: DP 5–6 at 1X, DP 6–9 at 2X and DP 9–10 at 3X. In contrast, GATK did not reach a plateau, as its accuracy continued to increase at higher thresholds, although the rate of improvement diminished as threshold rose (Fig.4.A and Additional file 10). Overall, the peak accuracy for Freebayes was observed at 3X DP 10, while GATK’s best performance was at 1X DP 10. On the contrary, the observed trend in average sensitivity was consistent for both SNP callers, showing a decrease at higher depth of coverage thresholds and lower sequencing coverage (Fig.4.B and Additional file 10). In terms of performance, Freebayes outperformed GATK in sensitivity (Additional file 10). At 1X, Freebayes average sensitivity varied from 9.88% at DP 1 to 2.85% at DP 5 and 0.54% at DP 10, whereas GATK sensitivity ranged from 8.50% at DP 1 to 1.83% at DP 5 and 0.51% at DP 10. At 5X coverage, Freebayes achieved an even better performance than GATK with 2.30% higher sensitivity values at DP 1 (35.27% vs. 32.97% with GATK), 4.30% at DP 5 (30.45% vs. 26.16%) and 4.22% at DP 10 (15.36% vs. 11.14%) (Fig.4.B). Among the common variants identified between the lcWGS datasets and the corresponding reference SNP panels, genotypic concordance was influenced by both SNP caller and sequencing parameters (Additional file 10). Freebayes consistently showed higher genotypic concordance under low coverage conditions compared to GATK (Fig.4C). For instance, at 1X coverage, Freebayes peaked at over 93.86% at DP 7, while GATK lagged behind, reaching a maximum concordance of around 85.90% at DP 6 for the same coverage. At higher coverages (5X), Freebayes maintained superior genotyping concordance, achieving a maximum of 98.72% at DP 10, while GATK’s concordance reached 1.45% less than Freebayes (Fig.4C). Discrepancies between true positives and the reference SNP panels genotypes were primarily due to heterozygous loci being misclassified as homozygous in the lc-dataset (Fig.4D). This error decreased with increasing minimum depth of coverage thresholds for both callers. With Freebayes, the percentage of misclassified heterozygous loci fell below 5% starting at DP 5 across all coverages, whereas GATK required higher coverage (3X and above) to achieve similar levels of concordance (Fig.4.D). Homozygous misclassifications were minimal for both callers (Fig.4E). Validation of LcWGS results with MAGIC lines As a proof-of-concept validation of the benchmark performed with the MEGGIC founders, the assessment was extended to four S5MEGGIC lines (Fig.1B), characterized by an intricate mosaic genome background from the MEGGIC founders (Additional file 1). Similarly, the four S5 lines were sequenced at 5X and downsampled at 4X to 1X with five replicates (Fig.1B). For the SNP calling, Freebayes was used, guided by insights obtained from the benchmarking analysis of the founders’ data. To assess the extent of potential false positives identified at each lc-DP combination, a comparative analysis was conducted using a gold standard as a reference (Fig.1B). The GS comprised a combination of shared and private founder variants from the 20X resequencing dataset, obtained by performing a Freebayes SNP calling at the population level (Fig.1C). This process yielded a total of 17,069,371 biallelic SNPs, considered reliable polymorphisms due to high read support, with their distribution provided in Additional file 11. Shared polymorphisms between S5MEGGIC lc-datasets and the gold standard (i.e., true positives) are the most informative in determining the most convenient lc-DP combination to maximize characterization accuracy and resource optimization. As expected, the number of true positives identified increased with the sequencing coverage and decreased at higher depth of coverage thresholds (Fig.5 and Additional file 12). Similarly to the founders’ benchmark (Fig.3), true positives decrements decelerated at higher sequencing coverage. Fixing DP 1 as 100% of TP for each sequencing coverage, 5X DP 10 still retained 21.67% of the TP versus 8.26% at 3X DP 10 and only 1.59% at 1X DP 10 (Additional file 12). Nevertheless, the proportion of true positives relative to the total biallelic SNPs identified for each lc-DP combination (%TP) did not follow the same trend and was different for each sequencing coverage (Additional file 13). The %TP slightly increased at higher sequencing coverage and the plateau shifted at higher depth of coverage thresholds when adding coverage. So that, at 1X, the highest %TP was 51.76% observed with DP 3, while it was 61.70% at 5X DP 10 (Additional file 13). Thus, considerations should be given to whether applying higher DP thresholds is advantageous. While this approach may reduce false positives and increase confidence in calling heterozygous loci, it could also result in a lower proportion of true positives relative to the total biallelic SNPs identified. Additionally, missing data for each lc-DP combination was assessed, as it is a variable that highly impacts downstream analysis (Fig.5 and Additional file 14). SNPs with 100% missing data after the filtering step were excluded from the total SNP count used to calculate the percentage of TP with 0%, 25% and 50% missing data. At low depth of coverage thresholds (DP ≤ 7), sequencing at 5X Page 16 of 16Baraja-Fonseca et al. BMC Plant Biology (2025) 25:1125 45. Dell’Acqua M, Gatti DM, Pea G, Cattonaro F, Coppens F, Magris G, et al. Genetic properties of the MAGIC maize population: a new platform for high definition QTL mapping in Zea Mays. Genome Biol. 2015;16:1–23. 46. Peterson GW, Dong Y, Horbach C, Fu YB. Genotyping-by-sequencing for plant genetic diversity analysis: a lab guide for SNP genotyping. Diversity. 2014;6:665–80. 47. Kumawat S, Raturi G, Dhiman P, Sudhakarn S, Rajora N, Thakral V, et al. Opportunity and challenges for whole-genome resequencing-based genotyping in plants. In: Sonah H, Goyal V, Shivaraj SM, Deshmukh RK, editors. Genotyping by sequencing for crop improvement. John Wiley & Sons, Ltd.; 2022. pp. 38–51. 48. Wragg D, Zhang W, Peterson S, Yerramilli M, Mellanby R, Schoenebeck JJ, et al. A cautionary Tale of low-pass sequencing and imputation with respect to haplotype accuracy. Genet Sel Evol. 2024;56:1–19. 49. Jiang Y, Jiang Y, Wang S, Zhang Q, Ding X. Optimal sequencing depth design for whole genome re-sequencing in pigs. BMC Bioinform. 2019;20:556. 50. Bhattarai G, Shi A, Mou B, Correll JC. Skim resequencing finely maps the downy mildew resistance loci RPF2 and RPF3 in spinach cultivars Whale and Lazio. Hortic Res. 2023;10:uhad076. 51. Sapkota S, Zou C, Ledbetter C, Underhill A, Sun Q, Gadoury D, et al. Discovery and genome-guided mapping of REN12 from Vitis amurensis, conferring strong, rapid resistance to grapevine powdery mildew. Hortic Res. 2023;10:uhad052. 52. Saripalli G, Adhikari L, Amos C, Kibriya A, Ahmed HI, Heuberger M, et al. Integration of genetic and genomics resources in Einkorn wheat enables precision mapping of important traits. Commun Biol. 2023;6:1–14. 53. Huang X, Feng Q, Qian Q, Zhao Q, Wang L, Wang A, et al. High-throughput genotyping by whole-genome resequencing. Genome Res. 2009;19:1068–76. 54. Liu J, Shen Q, Bao H. Comparison of seven SNP calling pipelines for the nextgeneration sequencing data of chickens. PLoS One. 2022;17:e0262574. 55. Martin AR, Atkinson EG, Chapman SB, Stevenson A, Stroud RE, Abebe T, et al. Low-coverage sequencing cost-effectively detects known and novel variation in underrepresented populations. Am J Hum Genet. 2021;108:656–68. 56. Watowich MM, Chiou KL, Graves B, Montague MJ, Brent LJN, Higham JP, et al. Best practices for genotype imputation from low-coverage sequencing data in natural populations. Mol Ecol Resour. 2023;00:1–13. 57. Barchi L, Rabanus-Wallace MT, Prohens J, Toppino L, Padmarasu S, Portis E, et al. Improved genome assembly and pan-genome provide key insights into eggplant domestication and breeding. Plant J. 2021;107:579–96. 58. Song K, Li L, Zhang G. Coverage recommendation for genotyping analysis of highly heterologous species using next-generation sequencing technology. Sci Rep. 2016;6:35736. 59. Kardos M, Waples RS. Low-coverage sequencing and Wahlund effect severely bias estimates of inbreeding, heterozygosity and effective population size in North American wolves. Mol Ecol. 2024;00:e17415. 60. Musich R, Cadle-Davidson L, Osier MV. Comparison of short-read sequence aligners indicates strengths and weaknesses for biologists to consider. Front Plant Sci. 2021;12:657240. 61. Schilbert HM, Rempel A, Pucker B. Comparison of read mapping and variant calling tools for the analysis of plant NGS data. Plants. 2020;9:439. 62. Wu X, Heffelfinger C, Zhao H, Dellaporta SL. Benchmarking variant identification tools for plant diversity discovery. BMC Genomics. 2019;20:701. 63. Barbitoff YA, Abasov R, Tvorogova VE, Glotov AS, Predeus AV. Systematic benchmark of state-of-the-art variant calling pipelines identifies major factors affecting accuracy of coding sequence variant discovery. BMC Genomics. 2022;23:1–17. 64. Stegemiller MR, Redden RR, Notter DR, Taylor T, Taylor JB, Cockett NE, et al. Using whole genome sequence to compare variant callers and breed differences of US sheep. Front Genet. 2023;13:1060882. 65. Ni G, Strom TM, Pausch H, Reimer C, Preisinger R, Simianer H, et al. Comparison among three variant callers and assessment of the accuracy of imputation from SNP array data to whole-genome sequence level in chicken. BMC Genomics. 2015;16:824. 66. Bhadhadhara K, Balamurugan M, Bharti N, Banerjee R, Kasibhatla SM, Joshi R. Performance Evaluation of Variant Calling Tools for Human and Microbial Genomes. In: 2023 International Conference on Emerging Trends in Networks and Computer Communications (ETNCC). Namibia: Institute of Electrical and Electronics Engineers Inc. 2023. pp. 235–42. 67. Bu M, Xu M, Tao S, Cui P, He B. Evaluation of different SNP analysis software and optimal mining process in tree species. Life. 2023;13:1069. 68. The 1000 Genomes Project Consortium. A map of human genome variation from population-scale sequencing. Nature. 2010;467:1061–73. 69. Espejo Valle-Inclan J, Besselink NJM, de Bruijn E, Cameron DL, Ebler J, Kutzera J, et al. A multi-platform reference for somatic structural variation detection. Cell Genomics. 2022;2:100139. 70. The 3.000 rice genomes project. The 3,000 rice genomes project. Gigascience. 2014;3:7. 71. 1001 Genomes Consortium. 1,135 genomes reveal the global pattern of polymorphism in Arabidopsis Thaliana. Cell. 2016;166:481–91. 72. Bukowski R, Guo X, Lu Y, Zou C, He B, Rong Z, et al. Construction of the thirdgeneration Zea Mays haplotype map. Gigascience. 2018;7:1–12. 73. Torkamaneh D, Laroche J, Valliyodan B, O’Donoughue L, Cober E, Rajcan I, et al. Soybean (Glycine max) haplotype map (GmHapMap): a universal resource for soybean translational and functional genomics. Plant Biotechnol J. 2021;19:324–34. 74. Jordan KW, Bradbury PJ, Miller ZR, Nyine M, He F, Fraser M, et al. Development of the wheat practical haplotype graph database as a resource for genotyping data storage and genotype imputation. G3 genes, genomes. Genet. 2022;12:jkab390. 75. Arrones A, Baraja-Fonseca V, Solana A, Plazas M, Soler S, Prohens J, et al. Resequencing and phenotyping of the first highly inbred eggplant multiparent population reveal SmLBD13 as a key gene associated with root morphology. Hortic Res. 2025;6:uhaf157. 76. Zhang W, Li W, Liu G, Gu L, Ye K, Zhang Y, et al. Evaluation for the effect of low-coverage sequencing on genomic selection in large yellow croaker. Aquaculture. 2021;534:736323. 77. Ye H, Ji C, Liu X, Bello SF, Guo L, Fang X, et al. Improvement of the accuracy of breeding value prediction for egg production traits in muscovy Duck using low-coverage whole-genome sequence data. Poult Sci. 2025;104:104812. Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.