BSA-seq was used to perform genomic mapping of virulence to three HF R genes: H6, Hdic, and H5. To do this, DNA bulks were prepared using DNA isolated from individuals segregating for virulence. This required that populations segregating for virulence and avirulence to each R gene had to be identified and that individuals within these populations had to be separately genotyped. Virulence to each wheat R gene examined presented different challenges. The solution was to prepare separate HF populations for each R gene. A description of each population is presented below, followed by a description of whole-genome sequencing and data analysis. Later, we describe the preparation of PCR-based markers that were used to improve genetic resolution.

The Hessian fly reference genome was used to map genes and identify PCR-based (microsatellite) markers and design PCR primers with the SSR Locator software (da Maia et al., 2008). The gene model identifiers are the names in the official HF gene set (OGS). These can be accessed at the USDA Arthropod i5k official workspace https://i5k.nal.usda.gov/data/Arthropoda/maydes-(Mayetiola_destructor)/GCA_000149185.1/ and in the genome assembly curated at the National Center for Biotechnology Information (NCBI), GenBank assembly accession number GCA_000149185.1. Molecular markers were used to genotype individuals taken from mapping populations and pooled DNA samples using standard PCR methods and the primers listed in Table S2 .

The end-sequences of HF genomic bacterial artificial chromosomes (BACs) that had been mapped to HF polytene chromosomes (Aggarwal et al., 2009) were used as part of the HF genome sequencing project (Zhao et al., 2015). We used these data to identify HF BACs that reside within HF genome scaffolds A1Random.66 and X2.8. From among these BACs, we selected BAC HF07L11 as a probe for scaffold X2.8 and BAC Md23L24 as a probe for scaffold A1R.66. Using methodology described previously (Stuart et al., 2014), these BAC clones were fluorescently labeled, denatured, and allowed to hybridize to complementary bases on HF polytene chromosome preparations. Later, the chromosomes were stained with DAPI and the positions of BAC hybridizations were examined and photographed using fluorescence microscopy.

To perform reverse transcription PCR (RT-PCR), total RNA was isolated from pools of 50 to 60 two-day-old first-instar larvae. RNA from avirulent and virulent strains were extracted separately using the RNeasy Mini Kit (Qiagen). Single-strand cDNA was reverse transcribed using the RNA of each individual pool separately and the SuperScript III First Strand kit (Invitrogen) according to the manufacturer’s recommendations. Single-strand cDNA was then used in PCR experiments using gene-specific primers ( Table S3 ). PCR was performed in 35 cycles of 95°C for 30 s, 55°C for 30 s, and 72°C for 60 s, followed with a final step of 72˚C for 5 min. The Hessian fly actin gene was included as a reference and internal control ( Table S3 ). RT-PCR products were visualized on agarose gels. RT-PCR amplifications were performed with at least three biological replications.

Phyre2 (http://www.sbg.bio.ic.ac.uk/~phyre2/) is a free web-based service for prediction of the three-dimensional (3D) structure of a protein sequence using homology modeling against a database of Hidden Markov Model (HMM) profiles of known 3D protein domain structures from the Protein Data Bank (PDB, http://www.wwpdb.org/) (Kelley et al., 2015). Phyre2 was used to examine the predicted protein structures of the two best candidate vHdic genes in an attempt to identify similarities with the domains of other proteins.

Total RNA was isolated from HF, first-instar, larval biotype GP using the RNeasy Mini Kit (Qiagen). RNA was subsequently reverse transcribed into first-strand cDNA using the SuperScript III First Strand Synthesis kit (Invitrogen). Double-stranded cDNA for effector genes was amplified from first-strand cDNA using gene-specific primers containing Gateway attB adapters ( Table S4 ). These primers were designed to exclude the corresponding secretion signal peptides from each effector gene. Gene attB PCR products were recombined into pENTR/pDONR vectors using the Gateway BP reaction (Invitrogen) and chemically transformed into Escherichia coli OmniMAX2 cells (Invitrogen). Recombinant colonies were selected on 50 µg/ml kanamycin LB plates. Colonies carrying the recombinant plasmids were selected for plasmid isolation and DNA insert sequencing. Genes in pENTR/pDONR vectors were recombined by Gateway LR reactions (Invitrogen) into expression vector pEDV6 and transformed into E. coli OmniMAX2 cells. Recombinant colonies were selected on 10 µg/ml gentamicin LB plates and used for plasmid isolation. A pEDV-GFP construct was built by LR recombination between plasmids pENTR1AGFP-N2 (FR1) (Campeau et al., 2009) and pEDV6 (Sohn et al., 2014). Plasmid pEDV5 (Sohn et al., 2007) was used as an empty vector (EV) control. pENTR1A-GFP-N2 (FR1) was a gift from Eric Campeau (Addgene plasmid 19364). Plasmids pEDV5/6 were generously supplied by J. Jones (The Sainsbury Laboratory, Norwich U.K.).

The pEDV constructs were transformed into B. glumae cells as follows: Electrocompetent B. glumae cells were prepared as previously described (Saitoh and Terauchi, 2013) with minor modifications. In brief, the B. glumae strain 336gr-1 was inoculated in 20 ml LB medium for 14-16 h at 28°C with shaking until OD600 = 0.8. The flask was then opened for 30 s under clean conditions and then maintained at 28°C with shaking for an additional 4 h. The cells were then pelleted twice at 4°C and 3000 rpm for 5 min and resuspended in cold 10% glycerol. The pellet was then dissolved in 600 μl of cold 10% glycerol and divided into 50-μl aliquots. These were stored at -80°C for transformations. Each plasmid construct (0.3 μg) was electroporated into B. glumae using a MicroPulse Electroporator (Bio-Rad). Transformant B. glumae strains were selected on gentamicin (25 µg/ml) LB agar. B. glumae strain 336gr-1 was a generous gift from J. H. Ham (Dept. Plant Pathology and Crop Physiology, Louisiana State University).

For HR assays, Bglu-pEDV strains were plated on King’s Broth (KB) agar 25 µg/ml gentamicin and incubated at 30°C for 14 to 16 h. Bglu-pEDV strains were dissolved in 0.9% NaCl solution at OD600 = 0.7. Bacteria suspensions were infiltrated with a 1-ml syringe without a needle into 4- to 5-week-old Nicotiana tabacum Burley 21 HA and N. benthamiana leaves. Infiltrated plants were maintained in a growth chamber at 24 ± 1°C with a 16:8 (light/dark) light cycle and 80 ± 5% relative humidity. HR was recorded after 24 h for N. tabacum and 48 h for N. benthamiana.

Bglu-pEDV cell suspensions were prepared and infiltrated in leaves of N. tabacum and N. benthamiana plants as described above. Leaf disks (150 mm diameter) were collected from the infiltrated areas using a cork borer 18-h post infiltration (hpi) for N. tabacum and 36 hpi for N. benthamiana. Leaf disks were floated on 15-ml nanopure water and incubated at 22°C with gentle shaking (100 rpm). Conductivity in the water was registered after 4 hours of incubation using a conductivity meter (Metler Toledo S30K) with a sensor probe (Conductivity Sensor LE703, Metler Toledo). Samples from three different plants were used as replicates for each treatment. Statistical analyses were performed using ANOVA and Tukey’s HSD for significant differences (p < 0.05) among treatments.

Zhao et al. (2015) mapped H6 virulence on the long arm of HF chromosome X2 using four different structured mapping populations and PCR-based DNA markers ( Table 1 ). This had been an arduous task. Therefore, we decided to test the capacity of BSA-seq in the HF using H6 virulence. H6-virulent bulk DNA was prepared from 23 males. H6-avirulent bulk was prepared from 19 males. Each male was collected from a field-derived population and genotyped in individual testcrosses ( Figure 1A ). These bulks were sequenced separately, resulting in 5.32 Gb of combined genomic data ( Table S1 ). The BWA tool (Li and Durbin, 2010) aligned the reads from each bulk against the HF reference genome and then SAMtools (Li et al., 2009) used these mapping data to identify 1.5 million SNPs. Popoolation2 (Kofler et al., 2011) calculated SNP allele frequencies within each bulk and then determined the differences in allele frequency for each SNP. Using these data, Fisher’s exact test (FET) was performed and average -Log10(FET) values within sliding 10-kb genome windows were plotted as a smoothed line across the HF genome ( Figure 2A ).

Unfortunately, the HF genome sequence is imperfectly assembled, and this was evident in the plotted data. Instead of a single peak that rose and fell over a single chromosomal position, several peaks were observed scattered along the genome map in both FET and FST plots ( Figure 2A, B ). These peaks indicate that linkage exists between H6 virulence and the underlying genome scaffolds. Most of these scaffolds are located in the “unmapped” fraction of the genome map. Their peaks suggest that they should be assigned to the chromosome known to carry the H6 virulence trait (Zhao et al., 2015), chromosome X2. Unmapped scaffold Un.18557, with an associated -Log10[FET] of 1.3, was the clearest example. Another peak with the same elevation was associated with a single chromosome A1 scaffold (A1R.66). The same logic suggests that all or part of scaffold A1R.66 also belongs to chromosome X2. This was later confirmed when Hdic virulence was mapped (described below).

Importantly, the peak with the greatest elevation (-Log10[FET] = 1.7) identified a 6.1-Mb region on the long arm of chromosome X2 ( Figure 2A ). Although this value failed to meet the statistical cutoff, established using the Bonferroni correction method (-Log10[FET] = 5.8), the 6.1-Mb genome region under this peak includes the 300-kb genome window where H6 virulence was previously mapped (Zhao et al., 2015). This region contains HF genome scaffolds X2.8, X2.10, X2.11, X2.12, and X2.13 ( Figure 2C ). Using Web Apollo (Kofler et al., 2011; Lee et al., 2013) and the HF genome reference sequence (https://i5k.nal.usda.gov/Mayetiola destructor) we identified 945 gene models within this window. Gene model Mdes009086-RA, which was previously identified as the candidate gene conditioning H6 virulence (vH6), resides within scaffold X2.11 ( Figure 2C ). Mdes009086-RA encodes a member of the putative HF effector protein family SSGP71. It is not transcribed in H6-virulent larvae (Zhao et al., 2015).

Although this investigation failed to resolve the position of H6 virulence as well as the previous investigation ( Table 1 ), the results demonstrated that non-structured field-derived populations can be used to identify SNPs linked to HF virulence and encouraged us to use BSA-seq to try to resolve the positions of virulence to other HF R genes. As discussed below, we believe that a larger mapping population would have improved H6 virulence mapping resolution.

To select a HF strain that was Hdic-virulent, the offspring of individual females taken from an Israeli field collection were examined for their ability to survive and stunt Hdic-wheat. During this process, we noted that all matings (n = 9) between individual Hdic-virulent females and individual Hdic-avirulent males produced Hdic-virulent male offspring. This suggested that Hdic virulence is X-linked because male offspring do not inherit X chromosomes from their fathers; matings between virulent females (v/v) and avirulent males (A/-) produce only virulent male offspring (v/-).

To test this possibility, we looked for linkage between X-linked microsatellite markers using conventional BSA. F₂ males collected from the Hdic virulence advanced intercross population (vHdic-AIP) were genotyped as Hdic-avirulent or Hdic-virulent individuals ( Figure 1B ). Separate avirulent and virulent DNA bulks were then prepared and used as template with primers designed for those microsatellites in separate PCR experiments. The two pools amplified alternative microsatellite alleles with six markers on scaffolds X2.7 and X2.8 (data not shown), indicating that those markers were linked to Hdic virulence. Genetic mapping using X2.7 and X2.8 markers and the DNA of 189 individual males collected from the F₂, F₆, and F₁₀ generations confirmed this linkage and suggested that Hdic virulence was present on either a 100-kb fragment on the proximal end of scaffold X2.8 or within the gap between X2.7 and X2.8.

For BSA-seq analysis, an Hdic-avirulent DNA bulk (n = 15) and an Hdic-virulent DNA bulk (n = 33) were prepared using DNA extracted from F₁₀ males ( Figure 1B ). Each bulk was separately sequenced. This produced 12.2 Gb of combined genomic data containing 1.2 million SNPs ( Table S1 ). The average FET values for SNPs across sliding 10-kb genome windows were plotted ( Figure 3A ). Two genomic regions were associated with peaks where -Log10(FET) values were greater than the statistical cutoff of 5.8. These appeared to be tightly linked to Hdic virulence. Consistent with the conventional BSA results, one region included a 1.1-Mb section of scaffold X2.8 ( Figure 3B ). Unexpectedly, the other region contained a 1.0-Mb section of scaffold A1R.66, the same A1 scaffold that showed X2 linkage in the H6 virulence investigation described above.

We hypothesized that A1R.66 was contiguous with scaffold X2.8, and present in the gap between scaffolds X2.8 and X2.7 ( Figure S1 ). That hypothesis was tested via genetic mapping using PCR-based A1R.66 markers ( Table S5 ) and 48 individual vHdic-AIP F₁₀ males. Consistent with the hypothesis, A1R.66 was linked to both X2 and Hdic virulence in those experiments ( Figure 3B ). It was tested again using fluorescence in situ hybridization (FISH) to the HF polytene chromosomes. An X2.8 BAC and an A1R.66 BAC used as probes hybridized together on the long arm of chromosome X2 ( Figure 3C ). Taken together, the genetic and FISH data determined that the most likely order of the scaffolds on the long arm of chromosome X2 is centromere-X2.7-A1R.66-X2.8. Therefore, disregarding the gap between scaffolds X2.8 and A1R.66, BSA-seq data positioned Hdic virulence within a contiguous 2.1-Mb region on the long arm of chromosome X2 ( Table 1 , Figure 3B ). The PCR markers used in this investigation further resolved the Hdic virulence, to within a 1.1-Mb region flanked by the closest recombinant markers, X2.8-202 and A1R66-169 ( Figure 3B and Table 1 ).

Previously, H5 virulence was mapped to a region spanning the centromere of chromosome A2. That region composed 30% of the chromosome’s length (Behura et al., 2004), and in three independent experiments, displayed severe recombination suppression. Therefore, in the present investigation, we first examined microsatellite markers identified on scaffolds flanking the A2 centromere for linkage to H5 virulence. To do this, we used the DNA extracted from 36 F₂ back-crossed males genotyped in the Behura et al. (2004) study ( Figure 1C ). The DNA of each male was examined separately for each marker. Consistent with the recombination suppression previously observed, each scaffold examined was linked to H5 virulence ( Table S7 ). However, recombination was completely lacking between H5 virulence and markers on scaffolds A2.6 (458.5 kb) and A2.7 (1.3 Mb), suggesting that the resolution of the position of H5 virulence might improve with a larger mapping population and additional markers.

In an attempt to do this, BSA-seq was applied to DNA bulks from 24 additional H5-avirulent and 24 additional H5-virulent (n = 24) males selected from the same F₂ back-crossed mapping population ( Figure 1C ). Because H5 virulence is autosomal, the H5-virulent bulk was developed from heterozygous (v/A) males while the H5-avirulent bulk was developed from homozygous (A/A) males. Therefore, SNPs in the H5-avirulent bulk were expected to be heterozygous in genomic regions unlinked to H5 virulence and homozygous in genomic regions linked to H5 virulence.

Whole-genome sequencing produced 13-Gb of high-quality reads, containing about 0.92 million SNPs ( Table S1 ). Average FET values for SNPs were plotted as described above ( Figure 4A ). Recombination suppression was again observed across much of the A2 chromosome as most A2 scaffolds had -Log10(FET) values of 2.0 and higher. In addition, one X1 scaffold and seven unmapped scaffolds had peaks with -Log10(FET) scores of 2.0 or greater, suggesting that they should probably be assigned to chromosome A2. Nevertheless, H5 virulence appeared to be most tightly linked to scaffold A2.7, with a statistically significant peak -Log10(FET) value of 8.0 ( Figure 4A ). Plotting -Log10(FET) values for 10-kb sliding-windows across the scaffold revealed that SNPs across the scaffold had similar FET values ( Figure 4B ). Consequently, we suspected that H5 virulence is probably associated with an Avr gene within the 1.3-Mb sequence of scaffold A2.7 ( Table 1 ).

Scaffold A2.7 contained 142 gene models. Thirty-seven of these encode predicted proteins carrying secretion signal peptides. BLASTP indicated that 26 of these are highly conserved and unlikely effector candidates. The remaining 11 genes have no similarities with genes in other insects. We discovered that one gene model, Mdes007142-RA, was composed of two putative effector-encoding genes ( Figure S3A ); one called SSGP47-1 (Chen et al., 2008) and another referred to here as Mdes007142(b). We examined the transcription of all 12 of these genes using RT-PCR on first instar larvae in the H5-avirulent and H5-virulent strains. None of these were differentially transcribed between bulks ( Figure S3B ). Therefore, multiple alignment and manual comparisons were performed on DNA sequences of each bulk against each other and the reference sequence. These comparisons examined the sequences of the virulence- and avirulence-associated alleles of each gene for evidence of frameshift and nonsense mutations that might disrupt the proper translation of virulence-associated sequences. These comparisons failed to identify putative null mutations that might be associated with the virulence allele of a candidate vH5 gene. Therefore, although BSA-seq was able to improve the resolution of H5 virulence on chromosome A2, we were unable to associate H5 virulence with a putative effector encoding gene. The possibility remains that a gene that has yet to be identified within the A2.7 sequence is responsible for H5 virulence.

Based on the similarities that HF-wheat interaction has with pathogen-plant systems, we decided to explore the pEDV system for bacterial T3SS-dependent delivery of effector proteins into plant cells and test whether HF candidate effectors are capable of suppressing plant defense responses. As described above, candidate vH6 (Mdes009086-RA) and one candidate vHdic (Mdes004160) were selected for this analysis. vH13 was included as a third HF Avr gene (Aggarwal et al., 2014). Green fluorescent protein (GFP) was used as a negative control. Each gene was moved separately into the pEDV6 vector ( Figure 5A ) and each construct transformed into B. glumae for infiltration in non-host N. tabacum and N. benthamiana plants (Kang et al., 2008; Sharma et al., 2013). At 24 hpi in N. tabacum and 48 hpi in N. benthamiana ( Figures 5B, D ), the HR that B. glumae normally induces in Nicotiana was evident on leaf tissue infiltrated with B. glumae carrying an empty vector (EV), GFP and candidate vHdic. However, HR was reduced or absent on leaf tissue infiltrated with B. glumae harboring candidate vH6 and vH13. These results were confirmed in a separate experiment in which ion leakage was used to assess plant cell integrity and membrane damage to the plant tissues (Rolny et al., 2011) ( Figures 5C, E ). Lack of HR suggests that, like many plant effectors, candidate-vH6 and vH13 proteins interfere with plant immunity.