A total of 33 lines of B. juncea spanning various genetic and geographic origins were evaluated for resistance against TuMV isolate UK 1 (Table 1). Two of these lines, TWBJ14 and TWBJ20, were taken forward and developed into two first backcross (BC₁) mapping populations shown to segregate for resistance to TuMV UK 1. Both TWBJ14 and TWBJ20 material were of the root morphotype and originated from China. The TuMV-susceptible B. juncea line 060DH17 line was used as the maternal parental line during initial crosses with TWBJ14 and TWBJ20 sources of TuMV resistance to produce the F₁ populations. This line is TuMV UK 1-susceptible, a doubled haploid and therefore homozygous across all genomic loci. Routine practices were employed for all plant crosses and for establishing an informative BC₁ population following phenotyping of relevant self and filial populations. The TWBJ14 and TWBJ20 BC₁ mapping populations were developed by crossing F₁ (060DH17 ♀ × TuMV-resistant plants ♂) UK 1-susceptible plants (♀) to TuMV UK 1-resistant S₁ plants (♂) that had been generated from the original resistant plants of each line used for initial crossing.

Resistance complementation testing was performed to elucidate whether resistance-associated genes were shared between TWBJ14 and TWBJ20 material through selective crossing and phenotyping of progeny with TuMV isolate GBR 6. Material used for complementation testing were F₂ populations generated by selfing the progeny of reciprocal crosses between TuMV-resistant TWBJ14 and TWBJ20 S₂ plants.

Brassica rapa ssp. perviridis cv. Tendergreen (Mustard-Spinach; TGM) is a highly susceptible host of TuMV and was used to maintain TuMV isolates through serial inoculations, and as a TuMV-susceptible experimental control to validate viral inocula during all testing. A plant line from the pathotyping system developed by Jenner and Walsh (1996), Brassica napus line R4 possessing TuRB01-based resistance to TuMV along with back-inoculation and reinoculation strategies were employed to confirm the authenticity and stability of TuMV isolates during phenotyping.

All experimental plant material was grown from seed within an insect-proof glasshouse maintained at 18 ± 2°C. Natural light levels were monitored and supplementary lights used in daytime when levels fell below 600 W/m². Seeds were planted in FP7 pots (7 cm square) containing Pot and Bedding M2 compost (Levington; medium grade sphagnum moss peat, pH 5.3–6.0). Plants were subsequently transferred to successively larger pots (FP9 and FP11; 9 and 11 cm square, respectively) containing John Innes No. 2 potting-on compost (Erin; loam-based compost with added fertilizer, pH ~6.5) when normal plant growth was hindered by available nutrition.

Four TuMV isolates, UK 1, v35Tunos +5570 A > G (vVIR24), CDN 1 and GBR 6 were employed within this study. They represent the most common TuMV pathotypes, 1, 3, 4 and 4, respectively (Jenner and Walsh, 1996; Jenner et al., 2000). TuMV UK 1 was obtained from swede in England (Tomlinson and Ward, 1978; Walsh, 1989) and is unable to overcome the TuMV resistance gene TuRB01 (Walsh et al., 1999). TuMV vVIR24 is an engineered mutant of UK 1 containing a single nucleotide mutation, conferring the ability to overcome TuRB01-based resistance (Jenner et al., 2000). TuMV CDN 1 was obtained from rutabaga in Canada (Walsh, 1989). It can overcome TuRB01 and at least three additional sources of resistance to UK 1 (Jenner et al., 2002b). Like CDN 1, TuMV GBR 6 is a pathotype 4 isolate. It has, however, been shown to overcome TuRB03-based TuMV resistance (Hughes et al., 2003); limited sources of resistance have been identified to this isolate.

As described by Walsh (1989), TuMV challenges were performed by mechanical inoculation of test plants at the 2–5 true leaf growth stage. TuMV-infected TGM leaf tissue was homogenized in cold inoculation buffer (K₂HPO₄, 10 g/L; Na₂SO₃, 1 g/L) using a sterile pestle and mortar. The resulting mixture was then abraded using sterile muslin cloth onto the adaxial surface of target plant leaves dusted with fine carborundum grit (~36 μm; ThermoFisher Scientific, United Kingdom). The abrasion and lysis of cells provided TuMV access to epidermal/parenchymatous cells in a controllable and uniform manner. The newest inoculated leaf on each plant was marked in order to discriminate between TuMV-inoculated and new, uninoculated leaves; this facilitated assessment of systemic viral spread. Uninfected TGM leaf tissue was used to mock-challenge experimental control plants.

Turnip mosaic virus-associated symptoms were scored weekly through visual inspection for four weeks post-inoculation, as described by Jenner and Walsh (1996). Qualitative visual phenotypes comprised of, “0,” resistance with no TuMV-associated symptoms, “R,” resistance where infection was limited to inoculated leaves with no systemic infection, “+,” systemic TuMV infection or +_N, systemic infection with necrosis (Figure 1). TuMV infection and titer were confirmed following visual assessment, using an indirect plate-trapped antigen ELISA (PTA-ELISA), as described by Walsh et al. (2002). Equal sized unchallenged and challenged leaves from test plants were harvested, macerated and diluted 1:1 with 0.05 M sodium carbonate coating buffer (Na₂CO₃, 1.6 g/L; NaHCO₃, 3 g/L). Test samples and appropriate TuMV-infected and uninfected control material were pipetted into predefined duplicate wells of a 96-well microtiter plate then incubated for 12 h at 4°C. Standard PTA-ELISA methods were followed with the first/detection antibody EMA67, a mouse monoclonal antibody, used at a concentration of 1/500. EMA67 has been shown capable of recognizing all tested isolates of TuMV by targeting a conserved TuMV coat protein epitope (Jenner et al., 1999). The second/visualization antibody, a goat anti-mouse IgG conjugated to alkaline phosphatase (Sigma-Aldrich, United Kingdom), was used at a concentration of 1/2,000. After all unbound antibodies were removed, plates were incubated at room temperature with p-nitrophenyl phosphate substrate solution (1 mg/ml in 0.1 M diethanolamine, pH 9.8). The reaction between the alkaline phosphatase conjugated to the detection antibody and this substrate produced a color change, which was used to indirectly quantify relative viral titer within samples.

Genotyping of mapping populations was performed using two methods. A total of 218 TWBJ14 BC₁ plants (52 resistant/0; 105 susceptible/+; 61 susceptible/+_N), 158 TWBJ20 BC₁ (22 resistant/0; 136/+_N) and the four parental plants used for the two initial F₁ crosses were genotyped using a 90 K single nucleotide polymorphism (SNP) array containing informative polymorphisms spanning all “A,” “B,” and “C” Brassica genomes (unpublished, SNP array under development and provided by the Parkin group, Agriculture and Agri-Food Canada; AAFC). This was complemented by genotyping by sequencing (GBS) on 102 TWBJ14 BC₁ (41 resistant/0; 31 susceptible/+; 30 susceptible/+_N) plant samples performed using the method of Poland et al. (2012). DNA was normalized to 20 ng/μl and 200 ng digested with PstI and MspI restriction enzymes at 37°C for 2 h. Adapters were then ligated to the digested DNA fragments using T4 DNA ligase for 2 h at 22°C. The reaction was inactivated, and samples pooled in two sets (96 samples in one set with those remaining processed in parallel alongside other samples, totaling 96). After pooling, the two libraries were amplified with a short extension time of 30 s and then purified using a QIAquick PCR Purification Kit (Qiagen). An Agilent Bioanalyzer was used to confirm fragment size and quality of the libraries. Libraries were then quantified using the Kapa library quantification kit (Roche) and sequenced on an Illumina HiScanSQ (100 bp paired-end reads). An existing pipeline was used to demultiplex sequencing reads. Trimmomatic (Bolger et al., 2014) was used to trim the adapters, and remove short reads and poor-quality data. Leading and trailing bases with quality below 15 and reads shorter than 55 bp were removed. The trimmed sequence reads were aligned to a B. juncea reference genome under development (IAP Parkin, AAFC, unpublished) using Bowtie2 (Langmead and Salzberg, 2012) with “-local” and “-sensitive” parameters used and “-score-min” of “L, 0, 0.8”. A custom Perl script was used to extract the best unique hits and the resulting BAM files were used for variant calling. UnifiedGenotyper with standard parameters from the Genome Analysis Toolkit (McKenna et al., 2010) was used to call polymorphic loci.

In addition to mapping resistance, a necrotic phenotype observed in the TWBJ14 BC₁ population was mapped using the 90K SNP array. The analysis was based solely on TuMV-susceptible plants, as it was not possible to determine whether resistant plants possessed the QTL associated with the necrotic phenotype induced by TuMV infection.

All QTL statistical analyses were performed using R v3.4.4. The add-on package “R/qtl” (Broman et al., 2003; Arends et al., 2010) was used to filter, model, and visualize genotypic data using inbuilt algorithms and according to established QTL mapping methods (Haley and Knott, 1992; Jones et al., 2009; Semagn et al., 2010; Xu et al., 2017). Prior to linkage group generation, a range of custom filters were successively applied to identify reliable marker subsets without the need for imputation or reliance on traditional measures of marker quality, such as sequencing read depth (Supplementary Tables 1 and 2). Pairwise recombination fraction (RF) values were calculated between all final marker subsets and a logarithm of the odds (LOD) score calculated for each pair based on a likelihood test of RF = 0.5. Linkage groups were inferred using appropriate RF and LOD score values adjusted to a stage where an appropriate number of linkage groups/chromosomes were generated. Ordering of markers within linkage groups was performed independently of reference genome position-associated physical loci. All possible orders within a successively increasing subset of markers were considered and applied where the number of obligate recombination events were reduced. After all markers had been ordered, linkage maps were finalized by calculating inter-marker distances based on RF using the Haldane mapping function, assuming no meiotic crossover interference. In addition to estimating inter-marker distances, the Haldane mapping algorithm was used to calculate an associated LOD score based on a test of how dependable calculated distances were. Using a selection of potential/assumed genotyping error rates, this LOD score was used to indicate the likely true error of genotyping present within final marker subsets.

Due to observed phenotypes being qualitative (i.e., distinct TuMV resistance against susceptibility), appropriate nonparametric algorithms were used for all QTL analyses. Multiple QTL modeling was first used to predict the most likely number of QTLs and any interaction between them using a stepwise forward/backward search algorithm, assuming a maximum of four distinct QTLs succeeded by permutation-based testing of the suggested QTL model (1,000 permutations, α < 0.05). Suitable Haley-Knott, multiple imputation-based, and expectation-maximization QTL interval mapping (IM) algorithms were then used to identify marker-phenotypic variance associations. Algorithm-specific statistical likelihood LOD scores indicating the probability of any marker being associated with the phenotype undergoing mapping were determined for all markers. To complement QTL modeling, cofactors were also applied through standard composite interval mapping (CIM) techniques to infer whether multiple QTLs with minor effects were likely present (Jansen, 1993; Zeng, 1994). Genome-wide LOD significance thresholds were determined for all QTL mapping algorithms through 10,000 permutations of randomized phenotypic data and α < 0.05. QTL CIs of a 1.5 LOD drop from the peak marker-associated LOD score were identified, and QTL intervals defined by extending this interval to adjacent flanking markers of each linkage map.

From the 33 B. juncea lines screened against TuMV isolate UK 1, eight were identified with some level of resistance (Table 1). Two S₂ populations produced by selfing both TuMV-resistant TWBJ14 S₁ and TWBJ20 S₁ plants exhibited uniform resistance when challenged with TuMV isolates UK 1, vVIR24 and CDN 1 (Table 2). This indicated that both plant lines possessed broad-spectrum resistance against isolates representing the most prevalent pathotypes (1, 3 and 4) of TuMV. Phenotyping of F₁ populations (060DH17 ♀ × TuMV-resistant plants ♂) for resistance to TuMV UK 1 suggested that both sources of resistance were recessive (Table 3). Subsequent phenotyping of both TWBJ14 and TWBJ20 F₂ and BC₁ populations further reinforced this; the observed patterns of segregation between TuMV resistance and susceptibility did not differ significantly according to chi-square (χ²) analyses from expected ratios based on a model of two recessive genes (TWBJ14 F₂ χ² = 0.094, TWBJ14 BC₁ χ² = 0.150, TWBJ20 BC₁ χ² = 2.38 < χ²_0.05 = 3.84; Table 4). TuMV-susceptible experimental control TGM plants challenged in all of these experiments were clearly infected with TuMV. Throughout all phenotyping experiments, there was also no indication that TuMV isolates behaved inconsistently when inoculated to the TuMV pathotype 1-resistant R4 plant line.

Inspection of 90K SNP array-derived genotypic data suggested raw data for both TWBJ14 and TWBJ20 BC₁ mapping populations was of high quality. A total of 11,732 marker loci were successfully genotyped with relatively few failing to be genotyped across less than 90% of plants from each of the BC₁ populations assessed (2.1 and 3.4%, respectively, for TWBJ14 and TWBJ20). The vast majority of markers provided data across more than 99% of the BC₁ samples tested (95.4 and 71.5% of TWBJ14 and TWBJ20, respectively). However, these marker loci were not strictly polymorphic and/or necessarily informative.

In contrast to the 90K SNP array data, all TWBJ14 GBS marker loci were identified based on a minimum of one sample being polymorphic to any other sample of the same BC₁ population assessed, so theoretically should have provided a higher abundance of informative and germplasm-specific data. Substantially more loci were immediately identified as potentially useful using this methodology and due to GBS not relying on pre-determined marker loci of interest. A total of 47,845 SNPs or insertion/deletion events (InDels) were identified. Irrespective of these having the potential of being informative, further inspection of GBS data suggested near-complete absence of any useful data. Of all GBS-derived marker loci, 85.1% failed to be genotyped across more than 10% of plants tested; there were no notable patterns linking any discrete sample(s) or marker loci to genotyping failure. Of the remaining genotyped loci, where more than 90% of BC₁ samples provided some level of data, sequencing read depth was highly variable. Due to the original parents of the TWBJ14 BC₁ population not being available for processing by GBS, prospective markers could not be filtered based on observations of appropriate allelic variation or hetero/homozygosity of absent parental genotypic data. Approximately half of the TWBJ14 BC₁ samples genotyped via the 90K SNP array (n = 218) were assessed using GBS (n = 102), further confounding issues with the efficacy of any recombination-based marker filtering strategies. These were limited due to inherently fewer recombination events being captured.

Successive filtering of marker loci proved highly effective when based on a combination of initial qualitative parameters (e.g., for 90K SNP array data, removal of markers where the doubled haploid parental line 060DH17 was genotyped as heterozygous) followed by iterative and optimizable measures, such as removing genotypic data generated across <75% of samples (Supplementary Tables 1 and 2). Applying this reductive approach to filter prospective markers, the 90K SNP array data generated a final subset of 1,663 and 3,158 markers, respectively from the TWBJ14 and TWBJ20 BC₁ populations. Recombination captured across each BC₁ population was identified as a limiting factor. The removal of redundant markers, where calculated RF was near-identical and which would consequently provide no additional linkage map resolution, reduced the final subset of markers to 1,064 for TWBJ14 and 966 markers for TWBJ20 (Supplementary Table 1). Assumed genotyping errors employed with the Haldane linkage mapping function suggested an impactful error rate of 0.00% across both these final subsets of 90K SNP array data.

Comprehensive exploration of filtering strategies based on sequencing read depth were unsuccessful in producing a TWBJ14 GBS marker subset capable of allowing linkage map production based exclusively on recombination. The use of imputation-based methods aimed at remediating the impact of genotyping failure also proved ineffective when implemented at any stage of marker filtering. Ultimately, the lack of genotyping data from the original TWBJ14 source of resistance and parents of the BC₁ population also proved extremely adverse. The identification of heterozygosity/recombination events across BC₁ plants was still assessable, based on observations of consistent genotyping of discrete alleles across BC₁ samples. Using this approach and similar reductive filtering as for the 90K SNP array data, a subset of robust markers was successfully identified from TWBJ14 GBS data (Supplementary Table 2). This final collection of markers was approximately 1% of the SNPs/InDels initially identified, totaling 481 unique loci. It is of note that this subset included informative markers which were otherwise discarded when implementing filtering strategies based on sequencing read depth. A small but potentially significant genotyping error rate of 0.55% was predicted when the final marker subset was assigned linkage groups and markers artificially ordered according to the B. juncea reference assembly.

Assignment of linkage groups and genetic map construction were performed during marker screening and employed as a method for the removal of unreliable genotypic data. The recombination fraction and associated LOD score calculated between approximately 7,100 filtered 90K SNP array markers from both B. juncea TWBJ14 and TWBJ20 BC₁ populations allowed assignment of markers into linkage groups based on a maximum RF of 0.15 and minimum LOD of 8 (Supplementary Table 1). This set of values retained the maximum number of informative marker loci possible, while remaining stringent. The high quality of 90K SNP array genotyping enabled marker loci to be assigned to 18 linkage groups, representing the 18 chromosomes of B. juncea, with few markers remaining unassigned. A single marker from TWBJ14 and 12 from TWBJ20 SNP array datasets were removed during linkage group assignment. Ordering 90K SNP array markers exclusively according to recombination fraction generated linkage maps with near-identical orders to those where markers were artificially ordered according to physical B. juncea reference genome loci (Figure 2). The only exception to this was within chromosome B04 where a partial chromosome inversion event was suggested for both TWBJ14 and TWBJ20 germplasm. Distribution of markers and overall marker density was consistent across final genetic linkage maps generated from filtered TWBJ14 and TWBJ20 BC₁ 90K SNP array marker data. The mean number of markers assigned to each chromosome for the TWBJ14 BC₁ population was 59.1 ± 3.7 (SEM), and 53.7 ± 2.7 (SEM) for TWBJ20 with no apparent gaps in genome-wide coverage (Figure 2).

A linkage map was also constructed based exclusively on pairwise recombination calculated between the most reliable 481 TWBJ14 BC₁ GBS markers remaining after excluding those presenting non-identical recombination. Compared to the physical order of markers on the B. juncea reference genome assembly, this linkage map could not be considered strictly reliable (Figure 3). Incorrect marker ordering occurred across many linkage groups, although, the majority of marker loci ordered by linkage were equivalent to orders predicted by physical reference assembly positions. Despite the suggestion from preliminary analyses of TWBJ14 GBS data that no viable data was present, a potentially informative linkage map was produced using considered screening of highly recalcitrant genotyping data.

A slight disparity was observed in TuMV susceptibility during resistance complementation testing of F₂ crosses depending upon whether TuMV-resistant TWBJ14 or TWBJ20 S₂ plants were paternal or maternal parents of two F₁ populations used for F₂ production. The F₂ population produced from germplasm, where the maternal F₁ parent was a TWBJ14 S₂ plant exhibited clear segregation for TuMV resistance with visual phenotypes correlating to PTA-ELISA data (Table 5). A contrasting less distinct segregation of TuMV resistance was recorded, where TWBJ20 material was the maternal F₁ parent.

Irrespective of F₁ parental maternity/paternity, segregation for resistance recorded for F₂ populations suggested that one TuMV resistance-associated locus was shared between TWBJ14 and TWBJ20. Where TWBJ14 germplasm represented the maternal F₁ parent, the resultant F₂ population segregated for TuMV resistance (41 resistant/0; 51 susceptible/+); this was determined by chi-squared (χ²) analysis to not be significantly different from the ratio of 7 resistant: 9 susceptible expected, where one resistant locus was shared (χ² = 0.025, value of p > 0.05) and differing significantly (χ² = 91.1, value of p < 0.001), where no resistance loci were in common (predicted segregation ratio of 31 resistant/0: 225 susceptible/+). Where TWBJ20 germplasm represented the maternal F₁ parent, respective χ² values of 3.49 (value of p > 0.05) and 118.6 (value of p < 0.001) were determined, confirming that one resistance locus was shared between the two B. juncea lines (Table 5).

In this study, we report the characterization of two new sources of broad-spectrum, recessive resistance effective against representative isolates of the most common pathotypes (1, 3 and 4) of TuMV. Both resistances involve two recessive genes. Resistance complementation testing and QTL analyses showed that the resistance-associated locus on chromosome A06 (retr04) is shared between TWBJ14 and TWBJ20, retr05 is specific to the former and retr06 is specific to the latter.

To date, at least 20 TuMV resistance-associated genes/loci have been mapped across all species of the Brassica genus (Li et al., 2019; Palukaitis and Kim, 2021) with significant efforts to characterize some (Walsh et al., 1999; Jenner et al., 2002a; Hughes et al., 2003; Qian et al., 2013; Nellist et al., 2014; Shopan et al., 2017). No sources of TuMV resistance across any Brassica species have been mapped to the “B” genome, and a single weak quantitative resistance to one TuMV isolate has been identified and mapped on the “C” genome (Walsh et al., 1999). Due to the evolutionary history of the Brassica genus, valuable traits of interest already identified in diploid Brassica genomes have the potential to be exploited in amphidiploid Brassica species through modern genetic editing techniques, or via interspecific crossing/resynthesis. This approach was used to great effect when introducing clubroot disease resistance from B. rapa into B. juncea material, wherein natural resistance is not evident (Hasan and Rahman, 2018) and to introduce turnip yellows virus resistance from Brassica oleracea and B. rapa into B. napus (Greer et al., 2021). Potential therefore exists for the generation of enhanced TuMV-resistant B. juncea cultivars using resistances already identified on the “A” genome of other Brassica species via resynthesis. Deployment of resistances in such a manner is, however, more difficult and time-consuming compared to the use of traditional MAS during natural introgression of any resistance-associated gene(s) from similar B. juncea material.

Few sources of TuMV resistance capable of being deployed directly through MAS have been identified within B. juncea germplasm (Nyalugwe et al., 2015; Shopan et al., 2017). The TWBJ14 and TWBJ20 B. juncea TuMV resistances identified and characterized here are consequently of high value. The vast majority of current resistances to TuMV in brassicas are also dominant and specific to discrete strains/pathotypes of TuMV. Such sources of genetic resistance to plant viruses remain valuable, but have potential to be rapidly overcome by less prevalent resistance-breaking viral mutants/recombinants (Stuthman et al., 2007). This can render such resistances ineffective if deployed irresponsibly owing to shifts in natural viral strain prevalence across a region, as is considered increasingly possible due to changes in pest and pathogen distribution predicted to occur with increasing worldwide temperatures (Bebber et al., 2013). The three recessive QTLs identified in this study appear to provide broad-spectrum resistance, as has also been shown for the retr01 recessive resistance in B. rapa (Walsh et al., 2002). If QTLs identified here are shown to continue to be effective against the most prevalent TuMV pathotypes, 1, 3 and 4, then germplasm and molecular markers explored here represent valuable breeding resources for developing enhanced B. juncea cultivars with broad-spectrum TuMV resistance.

Segregation ratios for TuMV resistance in F₁, F₂ and BC₁ TWBJ14 and TWBJ20 populations established that all sources of resistance investigated here were recessive. As retr04 was mapped to chromosome A06 and retr05 and retr06 were both mapped to chromosome A03, and the only other mapped recessive resistances to TuMV in brassicas have been mapped to chromosome A04 (retr01/retr02; Nellist et al., 2014) and chromosome A01 (retr03; Shopan et al., 2017), it is clear the resistances reported here are new. It is currently unknown how these new resistances to TuMV may vary in instances where single resistance-associated loci are present.

Sequence information used to produce the 90K SNP array markers that flanked the highly significant QTL for the necrotic phenotype on chromosome A06 was homologous to specific A06 loci of a reliable Brassica A genome reference assembly (Belser et al., 2018). This enabled comparison to published TuMV resistance genes. The dominant resistance gene TuRB01 and its likely ortholog TuRB01b were both mapped to chromosome A06 (Walsh et al., 1999; Lydiate et al., 2014). Further evidence was provided to indicate that these genes are similar or identical alleles and that B. napus may have acquired TuRB01 from the B. rapa gene pool during Brassica speciation (Walsh et al., 2002; Lydiate et al., 2014). Unpublished work undertaken by collaborators refined candidates of TuRB01/TuRB01b to a bacterial artificial chromosome contig. The full sequence of this contig aligned to chromosome A06 of the Belser et al. (2018) reference genome with 99.8% homology. This region was contained within the A06 necrotic phenotype QTL interval determined here. This could suggest the presence of TuRB01b in the plant lines investigated here, originating from a natural hybridization between a B. rapa plant possessing TuRB01b and B. nigra.

In terms of the TuMV resistances identified in this study, reliable identification of candidate genes for retr04, retr05, and retr06 was not possible with currently available data. This was due to the large number of genes suggested within QTL intervals and because robust gene annotation/validation has not been performed on the unpublished B. juncea reference assembly used in this study. Previously identified recessive TuMV resistance genes retr01, retr02 and retr03 have all been shown as variants of components of the eukaryotic translation initiation complex (Nellist et al., 2014; Shopan et al., 2017). Most recessive plant resistances to potyviruses are due to natural mutations in eukaryotic translation initiation factor 4E (eIF4E) genes (Robaglia and Caranta, 2006). As with retr01-based resistance in B. rapa, it is possible that eIF4E/eIF(iso)4E allelic variants (which TuMV cannot use) could be candidates for retr04, retr05, and retr06. Shopan et al. (2020) used sequence data of 95 established eukaryotic translation initiation factors (eIFs) from Arabidopsis thaliana to identify homologs across various Brassica species, including 190 putative B. juncea eIFs (BjueIFs). Using the positions of these BjueIFs on the B. juncea cv. Tumida assembly developed by Yang et al. (2016) and alignment of the retr04, retr05, and retr06 QTL-flanking SNP array markers, BjueIFs within all QTLs were identified (Supplementary Figure 1).

The only BjueIFs suggested in the chromosome A06 QTL associated with TWBJ20 were BjuA022587 and BjuA023685; these are both proposed by Shopan et al. (2020) to be homologs of eIF2ß. Of these two putative eIF2ß genes, BjuA023685 was also contained within the A06 QTL associated with TWBJ14 resistance. BjuA023685 may perhaps represents a candidate for retr04 shared by TWBJ14 and TWBJ20. Five additional BjueIFs (putative homologs of eIF2C, eIF3, eIF3δ, eIF4A, eIF4B, and eukaryotic release factor, RF) were also contained in the TWBJ14 A06 QTL (Supplementary Figure 1). Within the first A03 retr05/TWBJ14-associated QTL, six putative BjueIFs are suggested, all of which are subunits of eIF3, eIF4, or an eIF4A-like RNA helicase. Only two BjueIFs (BjuA010155; eIF4A and BjuA010508; eIF3F) were located within the retr06/TWBJ20 A03 QTL. While BjueIFs suggested in these QTLs may represent candidates for retr04, retr05, and retr06, the B. juncea cv. Tumida genome assembly employed by Shopan et al. (2020) may not accurately represent eIF variation in germplasm explored here. Future analyses of putative eIFs using the unpublished B. juncea reference assembly based on the specific 060DH17 line used in this study may elucidate if these eIFs are likely candidate resistance genes of interest. Molecular markers/polymorphic loci identified here could also be investigated against second backcross (BC₂) populations as part of work intending to refine candidate gene intervals, whilst validating prospective markers for commercial use. Additional characterization of the TWBJ14 and TWBJ20 resistances is required, but these new resistances present promise as a genetic resource for helping to mitigate TuMV.

For reasons comprehensively reviewed previously (Davey et al., 2011; Deschamps et al., 2012; Poland and Rife, 2012; Kim et al., 2016), GBS is considered to be an informative, reliable, and highly versatile genotyping option. In circumstances where an appropriate or reliable reference genome is unavailable, the massive high-throughput parallel sequencing associated with NGS-based techniques may represent one of few viable options for identifying genome-wide SNP/InDels in original germplasm. The use of GBS here facilitated the development of a reasonably robust linkage map using the most informative subset of 481 TWBJ14 BC₁-specific markers. The unfiltered data used to identify this subset presented significant issues, many of which are commonplace for GBS. Data produced by GBS tends to contain large quantities of missing genotyping calls compared to alternative methods, whilst also exhibiting inherent and often considerable error rates (Furuta et al., 2017). Missing GBS data can be attributed to presence/absence variation in germplasm under investigation, polymorphic restriction sites impacting the generation of reduced representation libraries, and/or differential methylation across DNA (Poland and Rife, 2012). Technical issues that contribute to missing data include the complexity of reduced genome representation libraries, number of unique sequence tags and overall sequence coverage. These issues may have occurred during employment of GBS here.

Sequencing read depth/coverage is the predominant measure used for rapid screening of NGS-based data for marker quality. As a case study, the processing of GBS data here represents a situation not previously reported. Highly recalcitrant genotypic data has been used to produce results comparable with those derived from near-optimal 90K SNP array data without any implementation of coverage-based screening of prospective markers. The consideration that this was achieved, where TWBJ14 90K SNP array data was generated from more than double the number of GBS BC₁ samples rendered this output increasingly striking. Read coverage cannot be considered an absolute indication of reliability nor measure for whether potential marker loci may ultimately prove informative. This was demonstrated here wherein “traditional” read depth-based screening of TWBJ14 GBS data proved unviable for producing a marker subset for effective linkage map production. Based on any filters applied according to sequencing coverage, a large quantity of GBS marker loci were considered reliable yet proved detrimental for effective linkage map production. The majority of these “reliable” marker loci presented qualitatively impossible genotyping because of the backcrossing strategy employed and allelic variation observed across BC₁ samples.

Other methods including the imputation-based methods of the “ABHgenotypeR” R package (Furuta et al., 2017) and hidden Markov models (Technow and Gerke, 2017) proved ineffective for processing the imperfect GBS data explored here. In contrast to using imputation-based methods, all questionable GBS marker loci identified here were removed without any amendments whatsoever made to original data. The methods implemented here limited prospective markers available for linkage map production and QTL analyses, although this was considered to not have significant impact due to the limited levels of recombination observed across the 102 TWBJ14 BC₁ plant samples assessed. Reductive marker screening approaches are speculated to be advantageous relative to strategies reliant on data amendment. A key advantage of this reductive approach relative to imputation is the removal of erroneous QTLs caused by false remediation of genotyping calls implemented because of cumulative genotyping errors or sparsity of markers.

It is debatable whether GBS is ever an appropriate method to employ for trait-mapping studies involving a large number of mapping population samples, due to both cost-efficiency and data handling considerations. Genotyping by sequencing routinely generates a large quantity of potentially useful markers. A level of diminishing returns can exist, however, for studies reliant on biparental mapping populations. The addition of markers where such a population has reached a stage of marker saturation with all recombination events captured can be considered detrimental, does not provide increased QTL interval resolution and may convolute analyses (Poland and Rife, 2012). This was likely the case for GBS performed here. Irrespective of the 47,848 marker loci initially suggested as perhaps informative prior to any screening for marker quality, captured recombination ultimately proved a limiting factor for whether markers were informative. During the final GBS marker loci filtering stage, wherein any loci demonstrating near-identical recombination fraction were removed, the reliable 1,057 marker subset was reduced to the representative 481 used for QTL analyses (Supplementary Table 2). The many benefits of GBS for QTL mapping, including no prior information of a target species genome being required and the ability to identify germplasm-specific markers, mean that such methods are highly valuable (Kim et al., 2016). This is regardless of whether recombination represents a bottleneck for identifying informative markers. Ultimately, the use of similar reductive screening methods to those employed here without read coverage being implemented as a predominant measure of marker quality may prove effective for other similar biparental QTL mapping studies.