Homoeologous chromosome exchange explains the creation of a QTL affecting soil-borne pathogen resistance in tobacco.

Crop plant partial resistance to plant pathogens controlled by quantitative trait loci (QTL) is desirable in cultivar development programmes because of its increased durability. Mechanisms underlying such resistance are difficult to study. We performed RNA-seq analyses for tobacco (Nicotiana tabacum) nearly isogenic lines (NILs) with and without favourable allele(s) at Phn7.1, a major QTL influencing partial resistance to the soil-borne pathogens Phytophthora nicotianae and Ralstonia solanacearum. Based upon combined analyses of transcriptome-based sequence variation and gene expression profiles, we concluded that allelic variability at the Phn7.1 locus was likely generated from homoeologous exchange, which led to deletion of low-expressing members of the SAR8.2 gene family and duplication of high-expressing SAR8.2 genes from a different subgenome of allotetraploid tobacco. The high expression of endogenous Phn7.1-associated SAR8.2 genes was correlated with observed resistance to P. nicotianae. Our findings suggest a role for genomic rearrangements in the generation of favourable genetic variability affecting resistance to pathogens in plants.

Introduction

Partial resistance to plant pathogens is of interest to plant breeders and to researchers studying the evolutionary dynamics of plant and pathogen populations. Although partial resistance can be highly variable and is of a lower level than that typically controlled by gene‐for‐gene interactions, it is considered an important component of plant improvement programmes because of its general non‐race specificity and increased durability (Cowger and Brown, 2019; Li et al., 2020; Poland et al., 2009). Hundreds of dominant genes providing complete resistance have been described (Kourelis and Van Der Hoorn, 2018). Comparatively, little is known about genetic components, or quantitative trait loci (QTL), controlling partial resistance (Li et al., 2020). This is likely due to the diversity of potential mechanisms that might be involved and also due to genetic complexities that make dissection of contributing factors more difficult (Cowger and Brown, 2019; Li et al., 2020). Increased knowledge of nucleotide variability controlling partial disease resistance in plants might help reduce crop loss from plant pathogens.

A classic example of partial disease resistance in crop plants is displayed by cultivated tobacco (Nicotiana tabacum L.) against black shank and bacterial wilt caused by the soil‐borne pathogens Phytophthora nicotianae and Ralstonia solanacearum, respectively. Prior quantitative trait locus (QTL) mapping experiments in this species resulted in the identification of multiple genomic regions with variable effects against black shank and bacterial wilt (Drake‐Stowe et al., 2017; Vontimita and Lewis, 2012; Xiao et al., 2013). A QTL now designated as Phn7.1 was identified on N. tabacum linkage group 7 (Nt7) and was observed to have the largest effect against black shank in all three experiments (Ma et al., 2019) and also against bacterial wilt in one investigation (Drake‐Stowe et al., 2017), making Phn7.1 a potential multiple disease resistance (MDR) locus (Wiesner‐Hanks and Nelson, 2016).

Because of the large effect of favourable allelic variability within the Phn7.1 region against multiple pathogens, there is interest in attempting to identify the nature of the DNA sequence variation responsible for this desirable effect. We previously transferred the Phn7.1 genomic region from highly resistant cultivar ‘Beinhart 1000’ to highly susceptible cultivar ‘Hicks’ using the backcross breeding method. The established BC4F3 nearly isogenic lines (NILs) and BC5F3 sub‐NILs are useful genetic tools because they isolate the allele(s) of interest in a common genetic background. These genetic materials possess as little as 12.6 cM of a Phn7.1‐associated region derived from Beinhart 1000 and exhibit significantly greater levels of field resistance to P. nicotianae (Ma et al., 2019). The objective of the research described here was to use RNA‐seq to (1) identify differences in gene expression between Hicks and a Hicks + Phn7.1 BC5F3 NIL with or without root exposure to P. nicotianae inoculum and (2) to determine if meaningful polymorphisms within genes of interest could be identified within the Phn7.1 region between these NILs. The overall objective was to identify potential candidate genes within the Phn7.1 region that are responsible for the favourable effect against disease progression caused by P. nicotianae.

Results

Identification of polymorphisms differentiating hicks and the hicks + Phn7.1 NIL

RNA‐seq data were generated for three to five biological replicates corresponding to each of the four different plant genotype × inoculation groups (GEO accession number GSE168854, and Table S1). Working towards the goal of identifying candidate genes underlying the favourable P. nicotianae resistance effect conditioned by Phn7.1, we first identified root transcriptome‐based SNPs and INDELs differentiating Hicks and the corresponding Hicks + Phn7.1 NIL (previously described BC5F3 NIL 235‐176 by Ma et al., 2019). The vast majority of identified polymorphisms were presumed to be associated with the introgressed Beinhart 1000 chromosome segment on Nt7, due to previous genome‐wide genotyping that determined the two lines to be nearly identical across the entire genome except for a chromosome segment associated with the Phn7.1 region on Nt7 (Ma et al., 2019). The vast majority of polymorphic sequences aligned to Nt14 of the current reference genome of tobacco cultivar K326, version Nitab4.5 (Edwards et al., 2017) (Figure 1a). This was surprising because three prior research studies mapped the major QTL affecting black shank resistance to Nt7 (Drake‐Stowe et al., 2017; Vontimitta and Lewis, 2012; Xiao et al., 2013). Microsatellite and SNP markers associated with this linkage group were used to create the Hicks + Phn7.1 NIL being studied (Ma et al., 2019). The majority of the identified potential SNPs appeared to be in the heterozygous condition in Hicks (an inbred line) and homozygous in the Hicks + Phn7.1 NIL (Figure 1b, Table S2).

Figure 1

Identification of Phn7.1‐associated region. (a) Chromosome map with SNPs identified for chromosome Nt7 (in blue circle) and chromosome Nt14 (in red circle). (b) Illustration of example SNP sites in modified snapshot of Integrative Genomic Viewer (IGV). Bar colours of brown, green, red and blue represent sequences of G, A, T and C, respectively. Boxed bars represent alternative sequences relative to the K326 reference genome. (c) Hypothetical composition of chromosomes Nt7 and Nt14 with respect to potential chromosome segment duplication due to homoeologous chromosome exchange between chromosomes of the S and T subgenomes, shaded in blue and red, respectively. (d) and (e) qPCR quantification of genomic sequences flanking two SNP sites in gene Nitab4.5_0000037g0020 and Nitab4.5_0000037g0170 for tobacco genotypes Hicks, Hicks + Phn7.1 NIL, and their F1 Hybrid. Red text and shading are associated with T subgenome (Nt14) genotypes. Blue text and shading are associated with S subgenome (Nt7) genotypes. Quantification scales are relative to Hicks (assigned a value of 1). (f) and (g) CAPS assays for SNP flanking sequences in tobacco genotypes indicated in qPCR assay above. Bands associated with S subgenome (Nt7) genotype are indicated with blue text, while bands associated with T subgenome (Nt14) genotype are indicated in red text.

We considered that the aforementioned observations might be explained by (1) the identified SNPs scored as heterozygous in Hicks are actually due to SNP variation between homoeologous loci of the different subgenomes of allotetraploid N. tabacum and (2) the Nt7 region in Hicks possessing the unfavourable allele(s) at Phn7.1 was replaced by the favourable resistance allele(s) ultimately derived from a homoeologous genomic region on Nt14 in P. nicotianae‐resistant lineages. N. tabacum chromosome Nt7 and Nt14 are reported to be homoeologous, originating from N. sylvestris (S genome) and N. tomentosiformis (T genome), respectively (Edwards et al., 2017). Consistent with this hypothesis, we found SNP alleles present for Hicks, but absent in the Hicks + Phn7.1 NIL and K326 (which carries the favourable allele at Phn7.1), to be present in N. sylvestris and absent for N. tomentosiformis (Sierro et al., 2013) and the current reference genome for N. tabacum cv. K326 (Figure 1b and Figure S1). We therefore considered it possible that a region carrying a P. nicotianae resistance allele(s) originating from the T genome became duplicated on a homoeologous chromosome segment from the S genome during evolution and/or breeding of some lineages of N. tabacum as shown in Figure 1c.

To gain further evidence for this hypothesis, we used quantitative PCR (qPCR) to study the relative copy number of genes within the region believed to have been duplicated or deleted in the Hicks + Phn7.1 NIL. These experiments were conducted for three genotypes: Hicks, the Hicks + Phn7.1 NIL, and a heterozygous F1 hybrid between the two. qPCR assays were conducted for SNPs at site Nt14:6750850‐2 in gene Nitab4.5_0000037g0020 and site Nt14:7355973‐5, which is the site of molecular marker Nt1AG1690 (Ma et al., 2019) in gene Nitab4.5_0000037g0170 (Figure 1b). These two sites allowed for the design of qPCR primers with at least two mismatches at the forward primer’s 3’ end for specific discriminatory amplification (Shi and Chiang, 2005) (Table S3). Consistent with the hypothesis of a duplicated/deleted genomic region (Figure 1c), the relative genomic DNA amounts for sequences of the S or T subgenome followed the patterns of 1 : 0.5 : 0 or 1 : 1.5 : 2 for the aforementioned three genotypes Hicks + Phn7.1 NIL, F1 hybrid, and Hicks, respectively. Some Nt7 sequences (originating from the S subgenome) do not exist in the Hicks + Phn7.1 NIL, and Phn7.1‐associated Nt14 T subgenome DNA sequences were present at a copy number of ~2× for the Hicks + Phn7.1 NIL compared to that for Hicks (Figure 1d, e).

Additional evidence of the potential for duplication/deletion was provided by CAPS marker assays for the sequences flanking these Phn7.1‐associated SNP sites (Figure 1b). Phn7.1‐associated Nt14 sequences (of T subgenome origin) can be detected in all tested tobacco lines with band density patterns congruent with observed qPCR patterns (Figure 1f, g). For example, the band density for Nt14 sequences was strongest for the Hicks + Phn7.1 NIL, and about twice that observed for Hicks. Nt7 (S subgenome) markers were not observed for the Hicks + Phn7.1 NIL. In contrast, genotyped S subgenome markers were observed for Hicks and the heterozygous F1 hybrid, with the marker band density being approximately two times greater for the former (Figure 2f, g).

Figure 2

Identification of candidate resistance genes within Phn7.1‐associated genomic region. (a) Illustration of Phn7.1‐associated region, including mapping densities for K326 short‐read sequences onto linkage group Nt14, distribution of Phn7.1‐associated scaffolds, previously reported Phn7.1‐associated SNP markers with respective significance tests reported as ‐log10(P) values (also known as LogWorth) in parentheses (Ma et al., 2019), SNPs identified in this study, and disease resistance‐related candidate genes in Phn7.1‐associated region. (b) Information for identified SAR8.2 genes from RNA‐seq reads, including assigned gene ID, subgenome assignment, partial sequence from first ATG and expression levels in roots illustrated using blue bars for Hicks and red bars for the Hicks + Phn7.1 NIL based on RNA‐seq‐based quantification in units of counts per million (cpm). (c) qPCR assays for copy ratio of S subgenome SAR8.2 genes (blue bars or text) and T subgenome SAR8.2 genes (red bars or text) in the genomic DNA of Hicks, Hicks + Phn7.1 NIL, and their F1 hybrid. qPCR based quantifications were relative to Hicks (set as 1). (d) CAPS assay for DNA specific to S subgenome SAR8.2 genes (bands indicated by blue text) and DNA specific to T subgenome SAR8.2 genes (bands indicated by red text) in the genomic DNA of Hicks, the Hicks + Phn7.1 NIL, and their F1 hybrid.

Further support for our hypothesis comes from investigation of DNA contents along the Nt14 linkage group of tobacco cultivar K326 that also carries the favourable allele(s) at Phn7.1 (Drake‐Stowe et al., 2017; Ma et al., 2019). Based on the mapping density of K326 short‐read sequences aligned to the current K326 reference genome (Nitab4.5), we found the DNA content of Nt14 in K326 to be apparently doubled in the region containing Phn7.1‐associated SAR8.2 genes, which is similar to the region defined by Phn7.1‐associated SNPs differentiating Hicks and the Hicks + Phn7.1 NIL (Figure 2a).

Identification of candidate resistance genes

An approximate ~6 Mb region of Nt14 of cultivar K326 defined by SNPs differentiating Hicks and the Hicks + Phn7.1 NIL contains six scaffolds collectively carrying greater than 100 putative gene models annotated in the genome assembly Nitab4.5 (Edwards et al., 2017) (Figure 2a ). Several additional genes were manually annotated. An additional group of 12 genes was found to be associated with polymorphic SNPs whose flanking regions were associated with two scaffolds unaffiliated with specific chromosomes in the current genome assembly (Table S4).

A number of these genes fell into one of three known groups associated with plant disease resistance. The first group contained several pathogenesis‐related (PR) genes. All of these appeared to reside outside of the genetic boundaries for Phn7.1 (Ma et al., 2019). A second group contained several tobacco mosaic virus (TMV) resistance‐like genes. These genes were not annotated in the current genome assembly. Some of these genes contained missense SNPs differentiating Hicks and the Hicks + Phn7.1 NIL (Table S4). There are no prior reports relating TMV resistance‐like genes to partial resistance against oomycetes or bacteria. A third group of genes were annotated as members of the SAR8.2 gene family, a small group of systemic acquired resistance (SAR) genes (Alexander et al., 1992; Shan and Goodwin, 2005). All four annotated SAR8.2 family members in the current K326 reference annotation could be identified within or around the Phn7.1‐associated genomic region. An additional SAR8.2 gene was manually annotated from a scaffold not assigned to a specific chromosome in the current reference genome assembly, but exhibits a Phn7.1‐associated SNP. All, or most, of the previously documented SAR8.2 gene family members appeared to be associated with the Phn7.1 genomic region (Figure 2a). Strong interest in the SAR8.2 gene family for potentially underlying observed Phn7.1‐mediated P. nicotianae resistance was supported by a previously issued patent claiming overexpression of members of this family for conveying increased resistance to P. nicotianae in N. tabacum (Ryals et al., 1998).

Investigation of the SAR8.2 gene family as it relates to Phn7.1

Members of the SAR8.2 gene family were considered as the mostly likely candidate genes contributing to increased P. nicotianae resistance in the Hicks + Phn7.1 NIL. Although four SAR8.2 genes are annotated in the current incomplete K326 reference genome (Edwards et al., 2017), the literature indicates that at least 15 SAR8.2 gene family members have been identified for N. tabacum (Alexander et al., 1992; Shan and Goodwin, 2005). By blasting previously known SAR8.2 transcript sequences against the K326 reference genome, more than ten loci containing SAR8.2‐like sequences were identified. These loci contained only incomplete or truncated SAR8.2 sequences. We therefore searched for all possible SAR8.2 transcripts within our RNA‐seq dataset. Fourteen distinct SAR8.2 transcripts, in addition to alternatively spliced versions, were identified from the combined root transcriptomes of Hicks and the Hicks + Phn7.1 NIL (Figure 2b and Table S5). Among these SAR8.2 gene family members, five were found to be specific to Hicks, with the majority featuring an ‘AT’ or ‘TT’ +17 and +18 nucleotides (nt) downstream of the ATG transcription site start. SAR8.2 gene family members with these ‘AT’ or ‘TT’ sequence tags are only present in the genome of N. sylvestris and are absent from the genome of N. tomentosiformis (Sierro et al., 2013). The remainder of the identified SAR8.2 genes are specific to the T subgenome, are present in Hicks and the Hicks + Phn7.1 NIL, and feature an ‘AC’ sequence tag +17 to +18 nt downstream of the ATG (Figure 2b).

The absence of Hicks or S subgenome‐specific SAR8.2 root transcripts in the Hicks + Phn7.1 NIL was consistent with the previously mentioned marker‐based findings suggesting the loss of a region of the Nt7 genome originating from N. sylvestris. On the other hand, the over 2× expression level of most T subgenome SAR8.2 gene family members in the roots of the Hicks + Phn7.1 NIL (Figure 2c) supports the hypothesis of a doubling in copy number of SAR8.2 gene family members originating from N. tomentosiformis in the genome of the Hicks + Phn7.1 NIL. qPCR results for the SAR8.2 gene family members in genomic DNA were also consistent with an over 2× increase in copy number for the N. tomentosiformis‐derived members in the Hicks + Phn7.1 NIL as compared to Hicks, and elimination of the N. sylvestris‐derived members in the Hicks + Phn7.1 NIL. SAR8.2 CAPS marker assays produced similar patterns. For example, marker bands corresponding to the T subgenome SAR8.2 genes were present for all studied tobacco lines, with band intensities being similar to patterns observed for the qPCR experiment. In contrast, marker bands for S subgenome SAR8.2 genes were only detected for Hicks and the Hicks + Phn7.1 NIL F1 hybrid (Figure 2d).

Associations between SAR8.2 genes and resistance to P. nicotianae infection

We used the amount of P. nicotianae RNA in inoculated samples (Figure 3a) as a measure of disease progression (Yan et al., 2017) to investigate the association between expression of SAR8.2 genes and resistance to disease progression. The amount of P. nicotianae in samples was represented by the overall alignment rate of RNA‐seq reads to the P. nicotianae transcriptome (Liu et al., 2016) for each sample (Figure 3a). To permit a comparison, SAR8.2 gene expression was also estimated using the overall alignment rate of RNA‐seq reads to all SAR8.2 transcripts for each sample (Figure 3b and Figure S3). The resulting expression profile patterns for all samples were confirmed by summing of RNA‐seq‐based quantifications (counts per million) for each SAR8.2 gene in the samples (Figure 3a,b, and Figure S3). Such overall expression analyses for SAR8.2 genes revealed the important observation that the overall level of SAR8.2 gene expression in control samples (uninoculated) of the Hicks + Phn7.1 NIL was found to be >2 times greater than that for the uninoculated Hicks control samples (Figure 3b,c, and Figure S3). This was primarily due to the high root expression of the SAR8.2 gene family members of the T subgenome, as compared to those of the S subgenome that accounted for only a very small portion of total SAR8.2 gene expression for Hicks (and Figure S3). According to the comparison, the overall expression of SAR8.2 genes in P. nicotianae‐inoculated samples was negatively correlated with pathogen transcript abundance (Figure 2a,b), inferring an additive/dosage inhibition effect against pathogen activity.

Figure 3

Functional association evaluation of Phn7.1 and SAR8.2 genes via RNA‐seq analysis. (a) Overall alignment rate of RNA‐seq libraries to the transcriptome of P. nicotianae and all identified SAR8.2 gene transcripts. (b) Expression of SAR8.2 genes based on RNA‐seq read counts in different RNA‐seq libraries mapped to each SAR8.2 gene transcript. All values were relative to the average of the sum of all SAR8.2 transcripts in Hicks (set as 1). The relative intensity of blue or red shades represents the expression level for S or T subgenome SAR8.2 genes, respectively. (c) Multidimensional scaling (MDS) plot for sample relationships based on the calculated variation distance (leading log fold change) between sample RNA‐seq libraries. Three distinct clustered groups were defined: control, resistant and susceptible (represented by black, blue and red circles, respectively). Refer to Figure 3a for sample IDs. (d) Upper panel of bar chart indicates the number of DEGs induced by Phn7.1 identified from a comparison between uninoculated Hicks and Hicks + Phn7.1 NIL control samples, and DEGs identified from the comparison between inoculated partially resistant or susceptible groups and all control group samples. Lower heatmap illustrates cross‐comparison of enriched GO terms for respective groups of DEGs illustrated in bar chart. Colour block in heatmap from yellow to red gradient represents different FDR levels; non‐colour in heatmap indicates non‐significant enrichment. Highly significant enriched GO terms (FDR < 1E‐20) are listed for down‐regulated DEGs specific for the ‘susceptible’ group vs control group comparison. (e) Relationships between overall alignment rates to P. nicotianae transcriptome and all SAR8.2 genes.

The relationship between SAR8.2 gene expression and P. nicotianae resistance was further evaluated in transcriptome‐based reactions of all samples to pathogen inoculation. Based on calculated expression variation distance analyses using EdgeR scripts (McCarthy et al., 2012) for RNA‐seq data generated for each plant genotype × inoculation combination, all samples could be clustered into one of three major groups (Figure 3c). One group included uninoculated (control) samples of Hicks and the Hicks + Phn7.1 NIL, where little total variation and fewer DEGs were observed amongst the samples (Figure 3c,d). The remaining two groups consisted of P. nicotianae‐inoculated samples. Of these, the ‘resistant’ group primarily consisted of the partially resistant genotype samples, Hicks + Phn7.1 NIL. The ‘susceptible’ group consisted primarily of the susceptible cultivar Hicks (Figure 3c). Due to the short time after pathogen inoculation, the susceptible group samples did not yet exhibit visible aboveground symptoms of disease. Disease‐associated signals were detected at the plant transcript level for these samples, however. For instance, gene functional analysis based on Gene Ontology (GO) enrichment for down‐regulated DEGs specific between the inoculated ‘susceptible’ and uninoculated control groups revealed a greater number of highly significant enriched GO terms related to protein synthesis, organelle activities, cytoskeleton structure, and many pathways related to basic cell activities (Figure 3d, and Table S6). Down‐regulation of genes associated with these GO terms could be expected to lead to reduced growth and even cell death, and is evidence of infection that was not yet visible between the ‘resistant’ group samples and uninoculated control samples. Accordingly, a higher overall alignment rate of RNA‐seq reads of the susceptible group samples to the P. nicotianae transcriptome was found as compared to resistant group samples (Figure 3d). It appeared that high overall expression of SAR8.2 genes was critical to the partial resistance response of the Hicks + Phn7.1 NIL in response to P. nicotianae infection (Figure 3e).

To gain further insight into the mechanism of partial resistance conferred by Phn7.1 or SAR8.2 genes, we investigated all 538 DEGs induced by Phn7.1 or high expression of SAR8.2s in uninoculated control samples of Hicks and the Hicks + Phn7.1 NIL. Only 10 DEGs from this analysis were located within the Phn7.1 genetic region (Figure 4a). These DEGs were all up‐regulated in the Hicks + Phn7.1 NIL, including three SAR8.2s which are highly expressed compared to other DEGs. Therefore, all other 528 DEGs, including those in unknown regions (within scaffolds not assigned to specific chromosomes) (Figure 4a) are more likely a consequence of the favourable genotype at the Phn7.1 locus, or due to the increased expression of SAR8.2 genes in the Hicks + Phn7.1 NIL.

Figure 4

Functional evaluation of Phn7.1 or SAR8.2‐induced DEGs in uninoculated control samples. (a) Genomic distribution of Phn7.1‐induced DEGs based on current K326 tobacco annotation (Nitab4.5). (b) Genes identified in biotic stress pathway via Mapman software for Phn7.1‐induced DEGs identified between uninoculated control samples of Hicks and the Hicks + Phn7.1 NIL. Boxes shaded in various intensities of blue and red represent up‐regulated and down‐regulated DEGs, respectively. (c) GO terms enriched for Phn7.1‐induced DEGs identified in control samples. ‘P’ and ‘F’ in the Ontology window represent biological process and molecular function, respectively. Sub‐terms are indented from primary terms in the Description section. The most highly indented terms are described at the gene function level. Lower FDRs are highlighted in red‐orange by increasing levels of intensity.

Through investigation of differentially expressed genes involved with biotic stress pathways using the Mapman platform (Thimm et al., 2004), we observed evidence of potential up‐regulation of pathogenesis‐related protein genes in the Hicks + Phn7.1 NIL vs Hicks, even in the absence of pathogen inoculation (Figure 4b, and Table S7). This finding is consistent with the results of overexpression of SAR8.2 transgenes in Arabidopsis, where AtPR‐1, 4 and 5 were up‐regulated (Lee and Hwang, 2006). This previous work did not examine, however, other PRs, such as PR‐10s (Fernandes et al., 2013), and PR‐3‐like endochitinase or glycoside hydrolase (Minic, 2008) that were identified as up‐regulated due to the high expression of endogenous SAR8.2s in this study. We also found other potentially disease resistance‐related genes to be up‐regulated in uninoculated control samples of Hicks + Phn7.1 NIL as compared to Hicks, including genes with leucine‐rich repeat domains (Padmanabhan et al., 2009), and genes encoding for S‐locus glycoproteins or like proteins (SGP/SGLP) (Figure 4c, and Table S6). In expanded GO enrichment analysis, most of the enriched GO terms for Phn7.1 up‐regulated DEGs in control samples overlapped with GO terms enriched for up‐regulated DEGs for inoculated ‘resistant’ samples (Figure 3d, and Table S6). Therefore, the existence of the favourable allele(s) at the Phn7.1 locus, which includes highly expressed SAR8.2 genes, can induce broader effects through up‐regulation of additional genes that might be involved with resistance to pathogens. Further strong evidence of the role of SAR8.2 genes on P. nicotianae resistance comes from field disease testing of tobacco lines expressing an RNAi construct designed to reduce expression of the SAR8.2 gene family, where resistance was observed to be very significantly reduced (Figure 5).

Figure 5

Increased susceptibility to P. nicotianae as the result of RNAi‐induced silencing of the SAR8.2 gene family in the genetic backgrounds of ‘Beinhart 1000’ and ‘K326’. (a) Histograms for area under disease progress curve (AUDPC) for RNAi lines and corresponding wild‐type genotypes. Means not sharing a common letter are significantly different from each other at the P < 0.05 level of significance. (b) Beinhart 1000 wild type and RNAi transgenic line P21a and (c) K326 wild type and K326 RNAi transgenic line 8a in P. nicotianae infested disease nursery. Photographs were taken 53 days after transplanting.

Discussion

Our findings indicate a relationship between increased copy number of highly expressed SAR8.2 genes and enhanced tobacco plant resistance to the soil‐borne pathogen P. nicotianae. Intergenomic recombination was likely responsible for creation of this genetic variation, whereby lowly expressed SAR8.2 gene family members located on Nt7 originating from N. sylvestris (S genome) were replaced by highly expressed SAR8.2 genes from a homoeologous chromosome of the T genome (ultimately derived from N. tomentosiformis). This structural variation led to a doubling of the copy number of highly expressed SAR8.2 genes in some tobacco lineages at the Phn7.1 locus.

Copy number variation (CNV) is an important type of genetic variation (Wellenreuther et al., 2019) that may have been previously underappreciated in populations. CNV has been shown to be associated with disease susceptibility in humans (Wellcome Trust Case Control Consortium, 2010) and many domestication traits in plants and animals (Lye and Purugganan, 2019). The mechanisms of CNV generation are not always clear, but transposable element activity is often considered as a major contributing force. For instance, evolution of copy number variation of multiple genes at the Rgh1 QTL locus has been associated with nematode resistance in soybeans (Cook et al., 2012) and is believed to be associated with past retrotransposon activity (Lee et al., 2015). In this report, we add to the body of knowledge describing homoeologous exchange as an alternative mechanism for evolution of CNV‐associated disease resistance QTL. Homoeologous exchange can occur in polyploid species to result in regions from one subgenome becoming duplicated/deleted to form presence/absence variation or segmental allopolyploidy (He et al., 2017; Hurgobin et al., 2018; Lloyd et al., 2018; Mason and Wendel, 2020; Stein et al., 2017). Quantitative disease resistance has previously been associated with presence/absence genetic variation (Gabur et al., 2018, 2020) and homoeologous recombination (Zhao et al., 2006) in Brassica napus. In addition, recent whole‐genome analyses using long‐read sequence data (Song et al., 2020) has revealed genetic variability affecting plant disease resistance to be enriched in genomic regions impacted by structural variation, including that due to homoeologous recombination in allopolyploid species (Chawla et al., 2021; Hurgobin et al., 2018; Samans et al., 2017).

SAR8.2 genes are one of several groups of genes whose expression is increased in response to pathogen infection, and that operate to trigger a heightened state of resistance against a broad range of pathogens in distal plant parts in a reaction known as systemic acquired resistance (SAR) (Alexander et al., 1993; Ward et al., 1991). To the best of the authors’ knowledge, the specific function of SAR8.2 gene products in mediating resistance remains unknown. The observed relationship between increased SAR8.2 gene expression due to a doubling in copy number of high‐expressing versions in tobacco genotypes carrying the favourable allele(s) at the Phn7.1 locus and disease resistance is consistent with reports describing similar effects in transgenic tobacco plants overexpressing a SAR8.2 gene (Ryals et al., 1998; Schultheiss et al., 2018) or transgenic Arabidopsis plants in which SAR8.2 genes from pepper were overexpressed (Lee and Hwang, 2006). Both overexpression of SAR8.2 transgenes in transgenic plants (Lee and Hwang, 2006) and increased SAR8.2 expression due to naturally occurring gene duplication was found to be associated with an up‐regulation of PR‐protein genes. Genome‐wide DEG identification and associated GO ontology enrichment revealed that high SAR8.2 expression was associated with up‐regulating the expression of additional groups of disease resistance‐related genes even in the absence of pathogen infection. Therefore, genes up‐regulated by favourable alleles at the Phn7.1 locus could provide additional layers of resistance to pathogen infection, in addition to that potentially provided by the SAR8.2 protein itself, which can exhibit antifungal and antibacterial activity (Lee and Hwang, 2006).

Due to the incomplete nature of publicly available genome assemblies for N. tabacum, many SAR8.2 genes were previously not well assembled and mapped to specific chromosomes. Based on our study, however, the majority of SAR8.2 genes appeared to be located within the Phn7.1 QTL region (Figure 2c,d). Additional investigation of SAR8.2 gene expression using existing RNA‐seq data from tobacco cultivar ‘TN 90’ (Sierro et al., 2013) also indicates that overall expression of SAR8.2 genes is largely root‐specific (Figure S3), which is consistent with observations that partial P. nicotianae resistance mediated by Phn7.1 is largely specific to roots and is not generally effective in aboveground tissues (McCorkle et al., 2013).

Experimental procedures

Genetic materials and inoculation with Phytophthora nicotianae

A BC5F3 NIL (designated as line 235‐176 by Ma et al., 2019) of black shank‐susceptible tobacco cultivar ‘Hicks’ was previously characterized as carrying an estimated 12.6 to 16.6 cM Beinhart 1000‐derived genomic region on Nt7 and exhibiting significantly greater levels of resistance to P. nicotianae as compared to Hicks (Ma et al., 2019). We first had the objective of carrying out RNA‐seq analyses on RNA extracted from inoculated and uninoculated roots of the genotype Hicks + Phn7.1 BC5F3 NIL (referred to hereafter as Hicks + Phn7.1 NIL) and Hicks to determine if informative differences in gene expression or meaningful polymorphisms could be identified between the two genotypes. Plants of each genotype were grown in sub‐irrigated 8.3 cm square pots filled with calcined clay (Turface Athletics All Sport, Profile Products, Buffalo Grove, IL, USA) and maintained at 27°C in a plant growth room with a 18 h light: 6 h dark photoperiod. The roots of each of eighteen 33‐day‐old plants of each genotype were inoculated according to McCorkle et al. (2013) with 10 mL of a 1 × 103 zoospores/mL suspension of a P. nicotianae race 0 isolate maintained at N.C. State University. A separate set of eighteen plants of each plant genotype was treated with 10 mL of the suspension buffer (without zoospores) and maintained as uninoculated control plants. All plants were maintained with sub‐irrigation for 42 h post‐inoculation, at which time roots were carefully rinsed to remove growth media, excised, transferred to liquid nitrogen for rapid freezing, and eventually stored at −80°C until RNA purification.

RNA‐seq data generation, processing and analysis

Frozen root samples of three to five biological replicates of each plant genotype × inoculation subgroup were ground into fine powder in liquid nitrogen using a mortar and pestle. Approximately 100 mg of ground root sample for each plant was transferred into 1.5 mL tubes for extraction of total RNA using a Qiagen Plant RNeasy mini kit (Qiagen, Hilden, Germany). RNA quality and quantity were checked using a NanoDrop UV‐Vis spectrophotometer phots instrument (Fisher Scientific, Waltham, MA) and Bioanalyzer (Agilent, Santa Clara, CA). Total RNA‐s were submitted to the Genomic Sciences Laboratory at North Carolina State University for construction of RNA‐seq libraries using the NEBNext® UltraTM Directional RNA Library Prep Kit for Illumina (New England BioLabs, Ipswich, MA) following the manufacturer’s protocol for long insertions. Prepared RNA‐seq libraries were pooled and provided to the Duke Center for Genomic and Computational Biology for sequencing on an Illumina HiSeq4000 instrument using the 150 bp paired end read mode. Reads were trimmed to remove poor quality bases using fastq‐mcf script with parameters of q = 30, l = 50 and P = 30. RNA‐seq libraries were mapped using HISAT2 (Kim et al., 2015) to the reference genome (version Nitab4.5) for N. tabacum cv. ‘K326’ (Edwards et al., 2017). Because some RNA samples were collected from inoculated plants, they were also expected to contain RNA from P. nicotianae. RNA‐seq reads were consequently also mapped to an available P. nicotianae transcriptome (Liu et al., 2016) using HISAT2.

Differences in levels of gene expression between root sample groups were investigated. Raw counts of RNA‐seq reads mapped to gene models were extracted using bedtools (Quinlan and Hall, 2010) based upon annotation for the N. tabacum cv. K326 genome annotation, version Nitab‐v4.5 (Edwards et al., 2017). Gene expression was analysed using the R package EdgeR (Robinson et al., 2010). Differentially expressed genes (DEGs) between sample groups were identified using classic options, with a false discovery rate (FDR) <0.05 and a fold change in expression >2. For specific gene family members of interest, respective expression levels were represented by numbers of RNA‐seq reads containing distinguishing nucleotide polymorphisms normalized for total read numbers for each RNA library. Expression levels in each sample were studied relative to the average of normalized reads of all three Hicks uninoculated control samples, which were set as 1 or 100% for efficient comparisons.

Transcriptome‐based DNA sequence polymorphisms distinguishing Hicks and the Hicks + Phn7.1 NIL were identified using a previously described pipeline (Shi et al., 2019). Sequence variation including SNPs or INDELs from RNA‐seq mapping files were extracted using the mpileup function of Samtools (version 1.6, Li et al., 2009). Phn7.1‐associated associated SNPs/INDELs differentiating Hicks and the corresponding NIL were further screened based upon the Genotype : Phre‐scaled Genotype Likelihood (GT : PL) section in the var file (version 4.2) generated by Samtools. To identify the most reliable SNP/INDELs, we selected those SNP/INDELs with the highest possibilities indicated by respective PL values to remain SNP/INDELs (Shi et al., 2019). For instance, SNP/INDEL sites GT : PL values of 0/0:0,255,255; 0/1:255,0,255; or 1/1:255,255,0 were selected as putative non‐sequence variation, heterozygous, or homozygous forms of sequence variation, respectively. All identified SNP sites were evaluated by visualization of respective RNA‐seq coverages aligned to the K326 reference genome for Hicks and the Hicks Phn7.1 NIL using Integrative Genomics Viewer (IGV, Robinson et al., 2011). In addition, IGV visualizations were also adopted for SNP sites in RNA‐seq read alignments of N. tomentosiformis and N. sylvestris (Sierro et al., 2013) for identification of variation origin (Figure 1b, and Figure S1).

Analysis of genomic DNA

Genomic DNAs were extracted from the genotypes of interest using filter paper‐based spin column method (Shi et al., 2018) and were used in Quantitative PCR (qPCR) and Cleaved Amplified Polymorphic Sequence (CAPS) analyses to determine the relative genomic copy ratio for genes thought be possibly associated with Phn7.1‐mediated partial resistance to P. nicotianae. The respective flanking sequence information for tested markers or SNPs within gene transcripts was extracted from RNA‐seq reads via ‘grep’ command in Linux. Primers for qPCR were designed to contain more than one mismatch with closely related sequences at the 3’ primer terminus to avoid non‐specific amplification (Shi and Chiang, 2005), and the gene elongation factor 1α (EF‐1alfa) (Schmidt and Delaney, 2010) was used as a reference for normalization of qPCR quantitation among different samples. CAPS Designer provided by ‘solgenomics.net’ website was used to design CAPS assays (Table S3) for genotyping of specific SNPs and identifying corresponding restriction enzymes. Mapping densities of illumine short‐read whole‐genome sequences of K326 (Edwards et al., 2017) to chromosome sequence in Nitab4.5 version genomic assembly was performed using HISAT2 in alignment, then using bedtools (Quinlan and Hall, 2010) to extract read counts in fixed‐size 100‐kb window along chromosome. Mapping densities were normalized using total reads of respective chromosome and calculated in unit of read counts per kb, which represents the roughly times of coverage based on 101 bp per read of DNA short‐read sequences.

Isolation and characterization of SAR8.2 gene family members

A family of SAR8.2 genes became of interest during the progression of the research described here. All published transcript sequences for N. tabacum SAR8.2 genes (Alexander et al., 1992; Shan and Goodwin, 2005) were used to query the N. tabacum cv. K326 reference genome (Edwards et al., 2017) and our transcriptome database to extract sequences with features common to the SAR8.2 gene family. Then two sequences, ‘GTAATATCCTCACAA’ and ‘TTTCTTTGGCTATTTTG,’ common to all reported tobacco SAR8.2s were used as ‘bait’ tag sequences in scripts to grep merged RNA‐seq paired end (PE) reads assembled by BBmerge script (Bushnell et al., 2017) of each RNA‐seq library. 268 nt‐length unique sequences from the ATG transcription start site were extracted and used for further extensions in additional rounds of grepping from merged RNA‐seq reads to obtain full‐length coding sequences. Naming of identified SAR8.2 genes followed published SAR8.2 sequences according to homology based on the CLUSTALW algorithm (Thompson et al., 1994). RNA‐seq‐based analysis of SAR8.2 gene expression was based on the number of RNA‐seq reads that aligned to the unique 268 nt‐length sequences identified via HISAT2 for the RNA‐seq libraries generated in this study, and published RNA‐seq datasets for N. tomentosiformis, N. sylvestris and N. tabacum cv. ‘TN 90’ (Sierro et al., 2013). Expression levels of the SAR8.2 genes in each sample were normalized by the total read number for the respective libraries. Overall expression of all SAR8.2 genes was also quantified by alignment rate to all identified SAR8.2 transcripts to avoid overestimation of RNA‐seq reads alignment to SAR8.2 genes with highly similar sequence.

Functional analyses

Suggested functional roles for identified DEGs in this research were determined using Mapman version 3.6.0RC1 (Thimm et al., 2004) and based on gene functional annotation prepared via the online website of the Mercator annotation tool version 3.6 using default parameters (Lohse et al., 2014). GO ontology enrichment analysis of DEGs was carried out using AgriGO version 2.0 (Tian et al., 2017) following the provided custom tools for singular enrichment analysis (SEA) and cross‐comparison of SEA (SEACOMPARE). GO enrichment was calculated using Fisher’s exact test with Yekutieli Hochberg’s multitest adjustment (FDR under dependency, P < 0.05) and the gene ontology type Plant GO slim.

RNA silencing of SAR8.2 gene family

Tobacco genotypes with Phn7.1, ‘Beinhart 1000’ (highly resistant to P. nicotianae) and ‘K326’ (medium resistance) were engineered with an RNAi construct designed to down‐regulate expression of the SAR8.2 gene family. The RNAi construct was prepared based on the pBI121 binary vector. A 115 bp segment of SAR8.2j exhibiting high sequence conservation with all SAR8.2 genes was selected as an RNAi silencing trigger sequence for simultaneous knockdown of all SAR8.2 genes (Miki et al., 2005). RNAi sense/antisense fragments and a 680 bp GUS linker fragment were amplified using primers with appropriate restriction sites at 3’ ends (Table S3), and assembled into sense : linker : antisense RNAi transgene fragments which were ultimately cloned into pBI121 by replacing the original GUS fragment.

Agrobacterium tumefaciens‐based transformation (An et al., 1986) was used to introduce the RNAi construct into the aforementioned two tobacco genotypes. Regenerated primary transformants (R0 individuals) were self‐pollinated in a greenhouse to produce R1 seed. Eight R1 individuals per R0 plant were grown in a greenhouse and testcrossed to their respective non‐engineered parental genotypes for the purpose of identifying R1 individuals homozygous for at least one transgene insertion. Non‐segregating testcross families were identified by germinating 100 seedlings on inorganic MS medium supplemented with 100 mg/L kanamycin. R2 seed was collected from plants determined to be homozygous for at least one transgene insertion.

Based upon a preliminary experiment, two non‐segregating, transgenic RNAi lines of the Beinhart 1000 background and three non‐segregating lines of the K326 genetic background were selected for evaluation of field resistance to P. nicotianae during the 2021 growing season. The five transgenic lines and their non‐engineered wild‐type parental lines were transplanted in a replicated experiment located in a field naturally infested with P. nicotianae and used for decades for evaluation of field resistance to this pathogen. The experimental design was a randomized complete block design with eight replications. Experimental units consisted of single 12‐plant rows. Within‐row and inter‐row plant spacing were 56 and 120 cm, respectively.

Initial stand counts were taken 26 days after transplanting. The number of plants killed by P. nicotianae was subsequently determined for each plot 22, 36, and 51 days after the initial stand counts. Percent death was determined for each plot at each disease rating, and the method of Madden et al. (2007) was used to calculate area under the disease progress curve (AUDPC). Data were analysed using PROC MIXED of SAS Version 9.4 (SAS Institute, Cary, NC), and mean separations were performed using Fisher’s protected LSD at α = 0.05 according to Steel et al. (1997).

Author contributions

R.S. and R.S.L. planned the research and designed the experiments. R.S. performed RNA‐seq experiments and analysed the data. J.J. and D.H.S. cultured the pathogen and performed inoculations. J.M.N collected field resistance data. R.S.L. and R.S. wrote the manuscript.

Conflict of interests

The authors declare no competing interests.

Figure S1 Snapshot of Integrative Genomic Viewer for SNPs identified on Nt7 and Nt14

Figure S2 Polygenetic trees of published SAR8.2s and SAR8.2 identified in the current study

Figure S3 Expression of SAR8.2s

Figure S4 Expression of Phn7.1‐induced DEGs in uninoculated Hicks and Hicks + Phn7.1 NIL samples

Table S1 RNA‐seq library information

Table S2 Phn7.1‐associated SNPs identified in current study

Table S3 Primer sequences

Table S4 Annotated genes in Phn7.1‐associated region defined by identified SNPs

Table S5 Information for SAR8.2 genes identified and studied in this investigation

Table S6 Results of Gene Ontology Analysis

Table S7 DEGs identified in biotic stress pathway via Mapman

Click here for additional data file.

Table S6 Results of Gene Ontology Analysis

Click here for additional data file.