The lines A25 and A4804 were used for the transcriptomic analysis. Line A25 is a selection of the market class fabada (white and very large seeds) that is susceptible to C. lindemuthianum race 38 (isolate Cl18). The NIL A4804 is a resistant genotype to C. lindemuthianum race 38 obtained from the cross A2806 × X4562. The NILs A2806 and X4562 are derived from A25 (Supplementary Figure 1), both having the seed phenotype of the market class fabada (Ferreira et al., 2012) and resistance to C. lindemuthianum race 38 controlled by a gene located in the Co-2 cluster (Ferreira et al., 2017). The resistance gene Co-2 is derived from Sanilac Bc6Are (Supplementary Figure 1), which was obtained from Cornell 49-242. Finally, sources of the Co-2 genes, Cornell 49-242 and Sanilac Bc6Are, were also included in this work.

The Cl18 monosporic isolate of C. lindemuthianum, classified as race 38 (Ferreira et al., 2008), was used in this work. To obtain abundant sporulation, the isolates were grown at 22°C in darkness for 10 days in potato-dextrose agar (Becton Dickinson and Company). Spore suspensions were prepared by flooding the plates with 5 ml of 0.01% Tween 20 in sterile distilled water and scraping the surface of the culture with a spatula. Inoculations were performed by spraying 10-d-old seedlings with a spore suspension containing 2 × 10⁶ spores mL⁻¹. Before sowing, seeds were disinfected in four steps: rinsed in distilled water to remove dirt particles, 30 s in 95% EtOH, 30 s in 15% hydrogen peroxide, and rinsed thoroughly in distilled water. The seedlings were maintained in a climate chamber at 23 to 24°C with 90 to 95% humidity and a 12-h photoperiod.

The experimental design had three replicates (corresponding to three resistance tests), two genotypes (susceptible genotype A25 and resistant genotype A4804), and three treatment assessment times: just before inoculation (named as 0), 24, and 48 h post-inoculation (hpi). On agar media at 24°C, the conidia germinated 4–6 h after sowing, and soon after, the appressoria were observed. At 48 h after sowing, there was extensive hyphal growth on the Petri plates, and a week later, sporulation was observed. Thus, the fungal attack started before 24 hpi, and the plant cells could detect the pathogen and start the cascade of reactions. In susceptible genotypes, after less than 24 hpi, the cytoplasm of infected cells gradually degenerates (O’Connell et al., 1985). Two seedlings per genotype were included in each replicate and treatment. The leaf tissues were harvested, flash-frozen in liquid nitrogen, and stored at −80°C before RNA extraction. In all, the study included 18 samples, named S (susceptible) or R (resistant), follow by the time when the leaf was collected (0, 24, or 48) and then the replicate number (1, 2 or 3). For example, S0.1 represents the susceptible genotype at 0 h/control from experiment 1.

Total RNA was isolated from samples using DNeasy Plant Mini Kit following the manufacturer’s instructions (Qiagen, Germany). RNA was quantified by fluorometric methods and quantity was investigated by using 2,100 Bioanalyzer Instrument (Agilent Technologies, United Kingdom). RNA libraries were prepared using the TruSeq Stranded mRNA Sample Preparation Kit (Illumina) and sequencing was carried out on the Illumina platform. The reads were mapped to the reference genome with HISAT2 splice-aware aligner (Kim et al., 2015) using the bean genome G19833 v1.0 (Schmutz et al., 2014).¹ Expression profiles are represented as read count and normalization values which were calculated based on transcript length and depth of coverage. The counts for mapped reads were normalized by calculating the FPKM (Fragments Per Kilobase of transcript per Million mapped reads). This analysis was performed in Macrogen Inc. (Seoul, Republic of Korea).

A principal component analysis (PCA) and hierarchical clustering analysis (HCA) were performed to detect the possible sources of noise in the results. The DEGs were identified through comparisons within resistant (A4804) and susceptible genotypes (A25) at 0 and 24 hpi and at 0 and 48 hpi (comparisons named as R24-R0, R48-R0, S24-S0, and S48-S0). In addition, the DEGs were investigated through comparisons between the two genotypes at 0, 24, and 48 hpi (named R0-S0, R24-S24, and R48-S48). The NOISeq package (2.38.0; Tarazona et al., 2011) and pheatmap 1.0.12 package in R project (R Core Team, 2021) were used to explore the quality of the samples and detect DEGs. The DEGs were identified using the criterion q > 0.80. Specific and common DEGs between genotypes and treatments were visualized using Venn diagrams constructed using the package ggVenn/ggplot2 in R project.

To investigate functional groups of DEGs in response to fungus, a GO analysis was performed using the Ensembl database (organism dataset: pvulgaris_eg_gene, version 2022-02-10) considering the three categories: biological process (BP), molecular function (MF), and cellular components (CC). The GO enrichment of significantly over-represented terms was carried out using the R package ViSEAGO (Brionne et al., 2019) and the DEG list. Enrichment tests were assessed using Fisher’s exact test (p ≤ 0.01) for the resistant and susceptible genotypes at both time level datasets [24 hpi (S0-S24, R0-R24) and 48 hpi (S0-S48 R0-R48)] and the three GO categories. The enriched GO terms were grouped on the basis of Wang’s semantic similarity into functional clusters using hierarchical clustering between GO terms with GO graph topology and Ward’s criterion.

The bean lines Cornell 49-242, Sanilac Bc6Are, A25, A2806, X4562, and A4804 were genotyped using the Genotyping by Sequencing method (Elshire et al., 2011) optimized in accordance with Schröder et al. (2016). DNAs were isolated from young leaves following the SILEX method (Vilanova et al., 2020), and DNA quality was checked in agarose gels. Genomic DNAs from the lines were digested individually with the Taqα1 and MseI restriction enzymes. Libraries were built based on the Poland et al. (2012) protocol with a different adaptor for ligation. Then, individual samples were checked by PCR, and the resulting products were visualized on 2% agarose gels. In total, 20 barcoded samples were pooled for PCR amplification. Sequencing was performed in the Illumina platform by Macrogen Inc. (Seoul, Republic of Korea). SNP calling was carried out by AllGenetics&Biology SL² using the P. vulgaris reference genome (genotypes G19833 v1.0,³ Schmutz et al., 2014). The SNPs supplied by GBS were filtered and extracted using TASSEL 5.1 software (Bradbury et al., 2007). SNPs with the following characteristics were considered in the analysis: (i) missing data, less than 50%; (ii) minor allele frequency, greater than 5%; (iii) mapped to one of the 11 pseudo-chromosomes. The genomic regions were delimited based on the physical position of the SNPs flanking the introgressed region (coordinates of the first and last introgressed SNPs). Finally, the GBS results reported by Murube et al. (2017) for the NILs A1258, A2806, and X2776 were used.

To identify variations useful to develop specific markers for the Co-2 region, raw RNA-seq reads of each experiment were assembled with the chromosome Pv11 of P. vulgaris v1 using the RBowtie package (Langmead et al., 2009). The alignments were visualized in Integrative Genomics Viewer software (Robinson et al., 2011), and the selected region was explored to detect polymorphism, insertions, or deletions (InDel), between the reference genome and the studied genotypes. The observed polymorphisms were also assessed by BLASTN and ClustalW alignments with the five bean genomes available (Phaseolus vulgaris G19833 v2.1; Phaseolus vulgaris 5 -593 v1.1; Phaseolus vulgaris UI111 v1.1; Phaseolus vulgaris Labor Ovalle v1.1).⁴ Those regions that showed an InDel in one genome were considered.

The primers for the PCR amplification of regions containing InDels were designed using the PrimerBlast tool (Ye et al., 2012). PCR reactions (20 μl) contained 2 μl of PCR buffer (10×), 1.2 μl MgCl2 (25 mm), 2 μl dNTP mixture, 0.8 μl of each primer (10 μm), 8.05 μl distilled water, 0.15 μl of TaKaRa LA Taq DNA Polymerase (5 U/μl; TaKaRa), and 50 μl 10 ng/μl DNA. The PCR reaction protocol performed on a Verity Thermal Cycler (Life Technologies Carlsbard, CA, United States) was as follows: an initial 5 min at 95°C; 30 cycles of 30 s at 95°C, 45 se at 62°C, 60 s at 72°C; final extension step at 72°C for 10 min. PCR products were separated on a 2% agarose gel, stained with RedSafe^™ Nucleic Acid Staining Solution (iNtRON, Seoul, Republic of Korea), and visualized under ultraviolet light.

At 7 days post-inoculation with C. lindemuthianum race 38, the susceptible genotype (A25) was dead, whereas symptoms were not observed on the resistant genotype A4804 (Figure 1). A total of 18 cDNA libraries were generated, one per genotype, replicate, and treatment, and, in total, they generated 876,756,658 clean reads using the NovaSeq platform (Supplementary Table 1). The reads of all the samples were used for transcriptome assembly, and an average of 91.1% of reads were mapped to the reference genome. Mapped reads were normalized by calculating the FPKM. A PCA of FPKMs estimated for each genotype, replicate, and treatment revealed two components that explained 79% of the variation. The biplot shows 3 samples, S0.1, S24.1, and R24.1, separated from the remaining 15 samples (Supplementary Figure 2A). Similarly, a hierarchical clustering analysis (HCA) classified the samples into two different main clusters separated from the samples S0.1, S24.1, and R24.1 (Supplementary Figure 2B). These samples had low RNA Integrity Numbers (RIN), less than 5.7, in the RNA extraction and were discarded from the analysis.

The DEGs were identified from seven comparisons (R24-R0, R48-R0, S24-S0, S48-S0, S0-R0, S24-R24, and S48-R48). A total of 5,740 differential expressions, involving 2,850 unique genes, were identified (Supplementary Table 2). A higher number of DEGs were found in the susceptible genotype than in the resistant genotype (Figure 2A) when the expression levels at 24 and 48 hpi were compared with the control (0 hpi). At 24 hpi after the susceptible genotype A25 was inoculated with C. lindemuthianum, 2,455 DEGs (686 upregulated; 1,769 downregulated) were detected, and the number decreased to 1,518 (469 upregulated; 1,049 downregulated) at 48 hpi (Figure 2A). Among the DEGs in the susceptible genotype, 1,615 genes only appeared in this genotype (Supplementary Figure 2B). The resistant genotype A4804 showed 831 DEGs (272 upregulated; 559 downregulated) at 24 hpi and 719 DEGs (306 upregulated; 413 downregulated) at 48 hpi. In addition, 127 DEGs were only detected in the resistant genotype, and 16 of them were identified at both 24 and 48 hpi (Figure 2B).

For the comparisons between resistant and susceptible genotypes at the same hpi, the highest number of DEGs, 143 (66 upregulated; 77 downregulated), were detected at 24 hpi. The number of DEGs decreased to 28 (16 upregulated; 12 downregulated) at 48 hpi (Figure 2C). In total, 21 genes were differentially expressed between both genotypes before inoculation and were discarded from the analysis. There were 14 DEGs between the resistant and susceptible genotypes that were common for the three hpi (Figure 2D). In contrast, there were 168 DEGs in both genotypes after inoculation, 10 of which maintained the differential expression at 24 and 48 hpi (Figure 2D): PHAVU_003G011800g, PHAVU_004G005400g, PHAVU_004G046400g, PHAVU_004G094000g, PHAVU_007G216700g, PHAVU_008G103500g, PHAVU_010G012900g, PHAVU_011G044000g, PHAVU_011G201700g, and PHAVU_011G203000g. Finally, four DEGs were only detected between both genotypes at 48 hpi: PHAVU_001G083000g, PHAVU_003G011000g, PHAVU_003G0939001g, and PHAVU_007G276500g.

A GO analysis was performed using the DEGs in each genotype. Among the 2,850 DEGs, 2,373 were assigned to at least one functional GO term in the Ensembl database (analysed 23/02/2022).⁵ To reveal the functional processes involved in the resistant and susceptible genotypes, the DEGs at 24 and 48 hpi were analyzed for enriched terms in the three GO categories (Supplementary Figures 3A–F, 4A–F). The differences were more evident at 48 than at 24 hpi, particularly in the BP category. The BP category contained 41 enriched GO terms at 24 hpi (Supplementary Table 4), 25 in the resistant genotype and 39 in the susceptible, with ‘translation’ and ‘biosynthetic process’ being the most significant pathways at both times. At 48 hpi, 111 GO terms were enriched (Supplementary Table 3), 89 in the resistant genotype, with ‘response to biotic stimulus,’ ‘cell wall organization,’ ‘hormone signalling,’ ‘flavonoid metabolism,’ ‘dephosphorylation,’ and ‘phosphatase activity’ being highly featured, and 31 in the susceptible genotype, with ‘translation’ and ‘biosynthetic process’ being the more prevalent terms. For MF at 24 hpi (Supplementary Table 4), only nine GO terms were enriched, four in the resistant and eight in the susceptible genotypes. In both genotypes, the term ‘structural molecule activity’ was the most enriched. At 48 hpi, 34 GO terms were significant (Supplementary Table 3) in MF, 31 in the resistant and 14 in the susceptible genotypes. Structural ‘molecule activity’ was the most enriched process in the susceptible genotype, whereas ‘hydrolase’ and ‘signaling receptor activity’ were the most enriched processes in the resistant genotype. For the CC category, 19 GO terms were enriched at 24 hpi (Supplementary Table 4), 13 of them in the resistant genotype and 14 in the susceptible genotype. The GO term ‘organelle ribosome’ was the most enriched in both cases, with ‘non-membrane-bounded organelle’ also being highly enriched in the resistant genotype. At 48 hpi 26 GO terms were enriched (Supplementary Table 3), 14 in the resistant genotype, with ‘extracellular’ and ‘cell wall regions’ being the most enriched term, and 15 in the susceptible, with ‘ribosome’ being the most enriched term.

Sequencing the GBS libraries yielded approximately 17.7 million reads in the six genotypes (Cornell 49-242, SanilacBc6Are, A25, A2806, X4562, and A2804), resulting in a total of 108,593 SNPs and 35,244 SNPs after filtering. The GBS revealed that 506 SNPs mapped on the end chromosome Pv11 (physical position>45 Mb). Genotypic comparisons of these 506 SNPs among the three resistant NILs harboring the Co-2 gene (A2806, X4562, and A2804), the resistance sources Cornell 49-242 and Sanilac Bc6Are, and the susceptible line A25 revealed that the three resistant NILs exhibited an introgressed region at the end of chromosome Pv11. These regions are tagged by SNPs having the donor resistance genotype (Sanilac Bc6Are). On the basis of the physical position of these SNPs, line A2806 has a region located at 45.10–48.78 Mb tagged by 183 SNPs of the SanilacBc6Are genotype (except for 21 isolated SNPs). The lines X4562 and A4804 maintained an introgressed region at 46.65–50.32 Mb that was tagged by 132 SNPs of the SanilacBc6Are genotype (Supplementary Figure 4). There was a common introgressed region among the three NILs located between 46.65 and 48.65 Mb.

The introgressed region in line A4804 carrying the Co-2 gene (46.65–48.65) has 165 annotated genes. Among them, 23 were differentially expressed in response to C. lindemuthianum. The representation of the differential expressions of the 23 genes in a heatmap revealed a separation between the resistant and susceptible genotypes (Figure 3) and the following four main groups of genes:

Group I, includes seven annotated genes that showed higher expression levels in the susceptible genotype and no changes in the resistant genotype, in response to C. lindemuthianum: *PHAVU_011G196600g, PHAVU_011G201800g, PHAVU_011G201700g, PHAVU_011G202000g, *PHAVU_011G193100g, PHAVU_011G203000g, and PHAVU_011G200600g. The genes marked with an asterisk encode hypothetical proteins with LRR domains.

Group II, contains five annotated genes that decreased in expression in response to C. lindemuthianum in both susceptible and resistant genotypes: PHAVU_011G200200g, PHAVU_011G192700g, PHAVU_011G200400g, PHAVU_011G200000g, and PHAVU_011G189900g.

Group III, includes eight genes with higher expression in the resistant genotype and non-changes in the susceptible genotype in response to C. lindemuthianum: PHAVU_011G193400g, PHAVU_011G193300g, PHAVU_011G194100g, PHAVU_011G193900g, PHAVU_011G193700g, *PHAVU_011G193800g, PHAVU_011G202200g, and *PHAVU_011G202100g. These genes have a serine/threonine-protein kinase function (6) or encode proteins with LRR domains (marked with an asterisk). They are located near two positions: a region with six genes (46.98–47.04 Mb, ‘core region’; Figure 4A) and another region with two genes (48.02 Mb).

Group IV, includes three annotated genes that tend to increase their expression in response to C. lindemuthianum, particularly in the susceptible genotype: PHAVU_011G206500g, PHAVU_011G2015000g, and PHAVU_011G2015001g. These genes have unknown functions in the bean reference genome.

The assembly of the raw transcriptome reads in the ‘core region’ in the chromosome Pv11 revealed InDels when read sequences of the resistant genotype (A4804) and the reference genome were aligned. The occurrence of InDel was also checked by BLASTN with the four available bean genomes and five of InDels were found in at least one of the genomes available. The susceptible genotype A25 did not show changes with the reference genome (G19833) in these four regions. The four polymorphic positions with InDel were:

M1 position, located in the third exon of the gene PHAVU_011G193000g (see Figure 4A). The resistant genotype A4804 has an insertion of 18 bp, also observed in the genomes 5–593 v1.1, UI111 v1.1, and, Labor Ovalle v1.1 (Supplementary Figure 5A).

M2 position, located in the intergenic region between the DEG PHAVU_011G193400g (Figure 4A) and PHAVU_011G193500g. The genotypes A4804, 5–593 v1.1, and UI111 v1.1 have a deletion of 10 bp. The reads with these sequences aligned with the gene Pv5-593.11G192900 annotated in the bean genome 5–593 (Supplementary Figure 5B).

M3 position, located in the gene PHAVU_011G193500g. The genotypes A4804 and 5–593 show a deletion of 6 bp and closed a mutation of two bp (Supplementary Figure 5C).

M4 position, located in the intergenic region between the genes PHAVU_011G193600g and PHAVU_011G193700g in the reference genome, but is annotated in the version 2 of G19833 in the gene Phvul.011G193600.1. The genotypes A4804, 5–593, and UI111 have a 12 bp deletion (Supplementary Figure 5D).

Specific primers were designed to amplify the four regions containing the InDels (Supplementary Table 5). PCR amplification generated fragments of the expected sizes, and the four markers were polymorphic between the resistant and susceptible genotypes A4804 and A25 (Figure 4B). The markers M1, M2, and M4 exhibited variations in size, whereas the M3 marker showed an amplification/non-amplification variation. These polymorphisms were also observed among the resistant genotypes Cornell49242, SanilacBc6Are, A2806, X2776, and X4562 (all of them with genotype as A4804), as well as the susceptible genotype A25.