The breeding line “ROR12” and the cultivar “Sprinter” were used in this study, together with a population of 148 F₇ RILs obtained from their F₂ progeny by single seed descent. “ROR12” was obtained from a local garden pea landrace by pure line selection (Pavan et al., 2016). Low strigolactone levels occurring in “ROR12” result in a slightly branched phenotype, which, however, does not cause a major penalty effect on the agronomic performance (Pavan et al., 2016). “Sprinter” is an old garden pea commercial cultivar previously shown to be highly susceptible to Oc (Bardaro et al., 2016; Pavan et al., 2016).

Two field trials (sowing dates 3 October 2020 and 8 January 2021) were carried out at the experimental farm “P. Martucci” of the University of Bari (41°01′22.1″N 16°54′21.0″E) in a silty–clayey experimental field continuously cultivated with legumes, known to be highly infested by Oc. RILs were arranged according to a randomized block design with three blocks and one replication per block, with each replicate consisting of 10 plants distant 0.15 m in a single row. The blocks were placed orthogonally to a gradient of Oc infestation observed in the previous 2 years. To check for the homogeneity of the Oc seedbank distribution within blocks, five replicates of the parental cultivar “Sprinter” were randomly allocated in each block as positive control. No fertilization and irrigation were applied during the growing season. Pest and pathogen management was performed using single applications of deltamethrin and difenoconazole, whereas weed control was performed with pendimethalin in pre-emergence and manual weeding in post-emergence. Genotypic response to Oc infection was evaluated at crop maturity on 25 May 2020, and 31 May 2021, as the average number of parasitic shoots emerged aboveground per plant.

Leaf tissue samples were collected from three individuals of the parental lines and one individual of each RIL. DNA was isolated using the DNeasy Plant Mini Kit (Qiagen) according to the manufacturer’s protocol and checked for quality and concentration using agarose gel (0.8%) electrophoresis and the Qubit 3.0 fluorometer (Life Technologies). A multiplexed ApeKI-GBS library was prepared as described by Elshire et al. (2011) and sequenced by a paired-end approach using the Illumina Novoseq 6000 sequencing system (Elshire Group Ltd.). After demultiplexing with the Axe algorithm (Murray and Borevitz, 2018), reads were trimmed for adapter and reverse-barcode sequences using the batch_trim.pl script from github (https://github.com/Lanilen/GBS-PreProcess). Alignment to the Pisum sativum v1.0 reference genome (Kreplak et al., 2019) was performed using bowtie2 (Langmead and Salzberg, 2012). After pooling together alignments of the biological replicates of the parental lines, the Stacks pipeline (Catchen et al., 2013) with the biparental filtering mode was used for SNP calling. Further filtering was performed in TASSEL 5.2.31 (Bradbury et al., 2007) by selecting SNP loci showing polymorphism between the parental lines and associated with call rate >90%, minor allele frequency (MAF) >0.25, and heterozygous calls<5%, and maintaining individuals displaying heterozygosity<10%.

The mstmap function of the ASMap R package (Taylor and Butler, 2017), implementing the Minimum Spanning Tree algorithm described by Wu et al. (2008), was used for linkage map construction. Default arguments within the function were applied, except for choosing the p-value threshold of 10^-11 for clustering marker loci in linkage groups. Marker loci associated with double crossover events or displaying skewed segregation (p< 0.05 after the Bonferroni correction) were removed from analysis using the drop.markers function of ASMap. The performance of the mstmap function was checked by plotting the heat map of pairwise recombination fractions (RFs) between markers and their pairwise logarithm of odds (LOD) score of linkage. Graphical representation of the linkage map was obtained using the iplotMap function of the qtlcharts R package (Broman, 2015). The relation between the linkage map and the pea reference genome was investigated by plotting genetic vs. physical positions for each chromosome, using the ggplot2 R package (Wickham, 2016).

Broad-sense heritability (H² _B) was estimated according to Schmidt et al. (2019), using the formula:

in which σ² _G is the genotypic variance, σ² _GY is the genotype-by-year variance, σ² _ϵ is the error variance, n_y is the number of years, and n_r is the number of replicates within each year. Restricted maximum likelihood (REML) estimates of variance components were obtained by fitting a random effect model with the lmer function of the lme4 R package (Bates et al., 2015), in which genotype, genotype-by-year, year, and replicates within years were set as random effects; the square root of phenotypic data was set as the dependent variable, with this transformation being necessary to correct for the right-skewness of the distribution. The model assumption of normality was assessed using the qqnorm function of the Stats R package (R Core Team, 2013). The ranova function of the lmerTest R package (Kuznetsova et al., 2017) was used to test the significance of the model random effect terms.

Mapping of QTLs was performed using linkage map data and best linear unbiased predictors (BLUPs) of RIL genotypic effects, which are widely used as alternative to phenotypic means to estimate genotypic values in QTL mapping studies (Ben Sadok et al., 2013; Allard et al., 2016; Molenaar et al., 2018; Wang et al., 2023). BLUPs were extracted by applying the ranef function of lme4 to the random effect model above described for H² _B estimation. The cim function of the Rqtl R package (Broman et al., 2003) was used to search for marker–trait associations by composite interval mapping (CIM). Arguments within the function were set to perform the Haley–Knott regression method and identify a LOD score QTL significance threshold based on a permutation test with 1,000 iterations. The plot function of the Stats R package (R Core Team, 2013) was applied on the output of the cim function to obtain graphs for chromosomal LOD scores. The percentage of the variance of RIL genotypic effects on the phenotype (PGE) explained by markers at QTL peaks was calculated with the formula:

Confidence intervals were identified using the lodint function of Rqtl, based on two LOD score units drop from the QTL peak. Genes included in the QTL confidence intervals were extracted from the pea reference genome annotation general feature format (gff) file, which was downloaded from the Unité de Recherche Génomique Info (URGI) bioinformatics platform (https://urgi.versailles.inra.fr/). The BLAST tools implemented by the same platform were used (a) to search, in QTL confidence intervals, for homologs of the strigolactone biosynthetic genes reviewed by Mashiguchi et al. (2021); (b) for physical mapping of the simple sequence repeat (SSR) marker AD174 (Loridon et al., 2005), previously associated with “ROR12” resistance (Bardaro et al., 2016); and (c) for physical mapping of the carbonic anhydrase gene Psat1g058960 that, according to the map of Carrillo et al. (2014), is closely linked to the RAPD marker OPAA19_702, in turn linked to the Oc resistance QTL n°br04 (Fondevilla et al., 2010). Search for defense response genes enrichment in QTL confidence intervals was performed by the gprofiler2 R package (Kolberg et al., 2020), using p = 0.05 as significance threshold and the Benjamini–Hochberg false discovery rate correction for multiple tests.

The effect of the different QTL combinations on the phenotype was investigated using the effectplot function of the Rqtl R package, which returned BLUP means and standard errors relative to different genotypic combinations at QTL pairs. Data were used to produce a custom plot using the ggplot2 R package (Wickham, 2016).

Kompetitive Allele Specific PCR (KASP) assays were performed using two allele-specific forward primers, marked with the FAM and HEX fluorescence dyes, and a common reverse primer ( Table 1 ). PCR reactions were performed at LGC genomics (Shanghai, China) according to standard protocols. Output fluorescence data were used to produce scatter plots, using the ggplot2 R package (Wickham, 2016).

Two-year trials were carried out in an experimental field known to be severely infested by Oc, aiming to evaluate the response of a segregant population of 148 F₇ RILs originating from a cross between “ROR12” and “Sprinter”. Scoring of the average number of parasitic shoots emerged aboveground per host plant indicated a fairly good uniformity of infestation within blocks, as relatively low dispersion around the mean was observed for five “Sprinter” replicates randomly allocated in each block (mean ± SD were 5.65 ± 0.35, 2.9 ± 0.14, and 2.95 ± 0.21 for the three blocks arranged in 2020, and 1.53 ± 0.12, 2.7 ± 0.46, and 3.4 ± 0.46 for the three blocks arranged in 2021).

The average number of parasitic shoots emerged aboveground per host plant ranged, for RILs, from 0 to 6.03 in 2020, and from 0 to 3.07 in 2021. In addition, in both years, this variable exhibited a distribution clearly deviating from normality, indicating the occurrence of one or a few major loci involved in resistance ( Figure 1 ). “ROR12” displayed very high resistance levels, with an average number of Oc shoots emerged per host plant of 0.13 and 0.06 in 2020 and 2021, respectively, suggesting the absence of transgressive segregation ( Figure 1 ). The estimated broad-sense heritability (H² _B) was 0.84, indicating a minor effect of environmental factors on phenotypic variation.

Sequencing of an ApeKI-GBS library obtained from the DNA of the parental genotypes and the RIL population resulted in approximately 2.2 million reads/sample. Approximately 33% of the reads were successfully mapped onto the pea reference genome (Kreplak et al., 2019). After the SNP call and quality control procedures, 6,182 polymorphic loci were identified. Further filtering to eliminate potentially spurious SNP calls, associated with loci displaying skewed segregation or associated with double crossover, resulted in a final panel of 4,489 markers. Of these, 4,127 were located on the seven pea chromosomes, whereas the remaining ones were located on pea superscaffolds and scaffolds. The number of polymorphic loci per chromosome and chromosome length displayed a moderate correlation (R ² = 0.55, p = 0.03), indicating a quite uniform distribution of variants across the genome.

Linkage analysis resulted in a genetic map containing seven linkage groups (LGs), in accordance with the pea haploid chromosome number ( Figure 2A ). Genetic length, physical length, and recombination rate associated with each linkage group are presented in Table 2 . Good collinearity was found between the position of markers in the genetic map and the one in the pea reference genome ( Figure 2B ). Most notable exceptions were represented by two regions of the chromosomes 4LG4 and 7LG7 ( Figure 2B ), and the mapping of 59 markers on a different chromosome than in the reference genome ( Supplementary Table S1 ). In addition, 362 SNP loci positioned on 9 superscaffolds and 125 scaffolds were anchored to the seven linkage groups ( Supplementary Table S2 ).

Testing for the effect of variance components indicated a non-significant contribution of the genotype-by-year interaction term ( Supplementary Table S3 ). Therefore, QTL mapping was performed using 2-year data, which were combined to obtain BLUPs of RIL genotypic effects on the phenotype (i.e., genotypic values). Three QTLs with LOD score peak above the significance threshold of 4.89 were identified. The QTL on chromosome 4LG4, named PsOcr-1, was associated with the highest LOD score peak (36.74), corresponding to 69.3% of the genotypic variance (σ² _G) ( Figure 3 ). The other two QTLs, located on chromosomes 1LG6 and 5LG3, were named PsOcr-2 and PsOcr-3, respectively. PsOcr-2 displayed a LOD score peak of 10.73 and explained 29.4% of σ² _G, whereas PsOcr-3 displayed a LOD score peak of 4.96 and explained 15% of σ² _G. The PsOcr-3 LOD score peak was mapped 14.98 Mb apart from the SSR marker AD174, previously associated with “ROR12” resistance (Bardaro et al., 2016). The PsOcr-2 LOD score peak was mapped approximately 17 Mb apart from the carbonic anhydrase gene Psat1g058960. This gene, according to the linkage map reported by Carrillo et al. (2014), is linked to the RAPD marker OPAA19_702, which was in turn linked to the Oc resistance QTL n°br04 by Fondevilla et al. (2010).

The QTL confidence intervals spanned physical regions of 17.2 Mb for PsOcr-1, 7.2 Mb for PsOcr-2, and 51.2 Mb for PsOcr-3, and contained 159, 99, and 617 genes, respectively ( Supplementary Table S4 ). Among these genes, 5 in PsOcr-1 and 25 in PsOcr-3 were annotated with the gene ontology (GO) term GO:0006952, “defense response biological process” ( Supplementary Table S4 ). Enrichment analysis showed significant enrichment of PsOcr-3 for this GO term (p-value = 7.08⁻¹⁰). The SNP loci corresponding to the PsOcr-1, PsOcr-2, and PsOcr-3 LOD score peaks were positioned within genes predicted to encode an F-box domain protein (Psat4g128720), a phenylalanine ammonia lyase (Psat1g046920), and a proton-dependent oligopeptide transporter (Psat5g221320), respectively ( Supplementary Table S4 ).

“ROR12” was previously shown to be a low-strigolactone line, which causes a reduced germination of Oc seeds (Pavan et al., 2016). Thus, we searched, within the QTL confidence intervals, for predicted genes showing homology with genes involved in the strigolactone biosynthetic pathway. This resulted in the identification, in PsOcr-2, of one 2-oxoglutarate-dependent dioxygenase (2OGD) (Psat1g046960) and, in PsOcr-3, of one cytochrome P450 oxygenase of the CYP711A subfamily (Psat5g201640), two cytochrome P450 oxygenases of the CYP722C subfamily (Psat5g209960 and Psat5g209880), and two 2OGDs (Psat5g206640 and Psat5g206800) ( Supplementary Table S4 ).

Information on the genotypic values and QTL genotypes of individual RILs was used to investigate the genetic effect of different QTL combinations. This indicated additivity between PsOcr-1 and PsOcr-2, and between PsOcr-1 and PsOcr-3 ( Figures 4A, B ), with the three QTL alleles contributing to Oc resistance all deriving from “ROR12”. Conversely, the genotype occurring at PsOcr-3, homozygous for either the “ROR12” allele (R) or the “Sprinter” allele (S), did not affect the genotypic value of RILs homozygous for the R allele at PsOcr-2 ( Figure 4C ), thus indicating epistasis of PsOcr-2 over PsOcr-3.

KASP assays were designed on the SNPs corresponding to the LOD score peaks of the QTLs PsOcr-1, PsOcr-2, and PsOcr-3. These were validated on the parental lines, as well as three different RIL panels, each one predicted from GBS data to include 10 lines homozygous for the “ROR12” allele, 10 lines homozygous for the “Sprinter” allele, and at least 1 heterozygous line. Each KASP assay yielded three fluorescence groups ( Figure 5 ). In addition, KASP genotypic calls were fully consistent with GBS calls.