Small RNAs generated by bidirectional transcription mediate silencing of RXLR effector genes in the oomycete Phytophthora sojae

Abstract

Oomycete pathogens secrete hundreds of effectors, including avirulence proteins that trigger host genotype-specific resistance response, to manipulate host immunity and facilitate infection. Sequence and expression variations of avirulence genes in pathogens are well known to be responsible for loss of host genotype-specific disease resistance. However, little is known on the underlying mechanisms associated with virulence variation in the diploid Phytophthora pathogens. We report in this study that the endogenous small RNAs (sRNAs) are involved in the variation of expression of avirulence gene Avr1b in P. sojae. The sRNAs were originated from the natural antisense transcripts of Avr1b. We further showed that the sense and antisense expression of Avr1b were programmed by the 10-base deletions in their promoter regions. Expanded analysis showed that up to 31% of the P. sojae RXLR effector genes were associated with sRNAs. Genome analysis further showed that the 9-bp and 10-bp insertion/deletion variants were significantly enriched in the regulatory regions of RXLR effector genes. These results indicate that the expression of RXLR effector genes are programmed by significantly enriched variations in their regulatory regions that lead to the variations in bidirectional transcription, which likely further affect production of endogenous sRNAs and silencing of homologous RXLR effector genes of Phytophthora pathogens.

Background

Oomycetes are eukaryotic pathogens causing severe damages to crops (Kamoun et al. 2015). For example, the worldwide annual losses caused by the Irish potato famine pathogen Phytophthora infestans is about $6.7 billion (Haverkort et al. 2008; Haas et al. 2009). The losses caused by soybean root and stem rot pathogen Phytophthora sojae is estimated at $1–2 billion per year (Wrather and Koenning 2006; Tyler 2007). Oomycete pathogens secrete a large array of effector proteins to manipulate host immunity and facilitate infection (Kamoun 2006; Stassen and Van den Ackerveken 2011). Many effectors, such as the RXLR family effectors (Jiang et al. 2008; Haas et al. 2009), are able to trigger or suppress programmed cell death in planta (Oh et al. 2009; Wang et al. 2011). If an effector is recognized by a specialized host resistance (R) protein, typically a nucleotide-binding site and leucine-rich repeat protein, it triggers hypersensitive response (HR) and prevents proliferation of the pathogen, and is called avirulence (Avr) protein (Jones and Dangl 2006). Interestingly, nearly all the avirulence genes identified in oomycetes encode RXLR effectors (Jiang and Tyler 2012; Anderson et al. 2015).

The successful pathogens have evolved diverse genetic strategies to evade the recognition by host plants. A variety of genetic variations in Avr genes, such as nucleotide substitution, insertion, deletion, frameshift, stop-codon mutation, truncation and copy number variation, are shown involved in the generation of novel types of virulent strains (Shan et al. 2004; Tyler et al. 2006; Qutob et al. 2009). Reversible virulence/avirulence variations of Phytophthora pathogens were frequently observed during asexual reproduction (Rutherford et al. 1985; Malvick and Percich 1998; Abu-El Samen et al. 2003). One possible reason is the suppressed expression of avirulence genes, thus leading to loss of race-specific resistance. P. sojae avirulence gene Avr1b can trigger HR in soybean plants carrying resistance gene Rps1b (Shan et al. 2004). Interestingly, epigenetic variation is also observed in the regulation of Avr1b phenotype in P. sojae. Two P. sojae strains, P6594 and P6497, carry the identical Avr1b gene, but displayed distinct phenotypes on Rps1b soybean plants. The accumulation of Avr1b mRNA in strains P6594 and P6497 are closely related to their Avr1b⁺ and Avr1b⁻ phenotypes, respectively. Avr1b of P. sojae strain P7064 (Avr1b⁺) has a synonymous single nucleotide polymorphism and is an epiallele of Avr1b^P6497 in fact. Genetic cross of strain P6497 and strain P7064 defined that the expression of P. sojae Avr1b is regulated by a trans-acting factor, and these two genes determine the Avr1b phenotype (Shan et al. 2004). However, how the mRNA level of Avr1b is regulated remains unclear.

Two popular hypotheses could address the epigenetic phenomenon at Avr1b locus: one is DNA methylation leading to transcriptional gene silencing, and the other is RNA silencing at posttranscriptional level and guided by small RNA (sRNA). However, it is known in oomycetes that transgene-triggered transcriptional gene silencing is not associated with DNA methylation (van West et al. 1999; Judelson and Tani 2007; van West et al. 2008). Two important gene silencing components involved in DNA methylation, cytosine methyltransferases and RNA polymerase IV, are also absent in the P. infestans genome (Vetukuri et al. 2011a). Moreover, many sRNAs are discovered in oomycetes and some of them were shown to be associated with silencing of RXLR effector genes (Vetukuri et al. 2012; Fahlgren et al. 2013; Qutob et al. 2013; Jia et al. 2017).

In this paper, we show that the silencing of Avr1b in P. sojae is associated with sRNA derived from its natural antisense transcripts (NATs). Close examination showed that the bidirectional transcription of Avr1b locus is likely programmed by the insertion/deletion (INDEL) variations in its regulatory regions. Genome analysis further showed that the 9-bp and 10-bp INDEL variants were significantly enriched in the regulatory regions of RXLR effector genes. Finally, analysis of the gene expression showed that up to 31% of RXLR effector genes in P. sojae were associated with sRNAs. These results indicated the expression of RXLR effector genes is likely programmed by the significantly enriched variations in their regulatory regions, leading to the variations in bidirectional transcription, production of endogenous sRNAs, and efficient silencing of homologous RXLR effector genes of Phytophthora pathogens.

Results

Identification of Avr1b-associated sRNAs

Avr1b of P. sojae is the first avirulence gene cloned from oomycete pathogens and is recognized by soybeans carrying Rps1b (Shan et al. 2004). Two types of virulent P. sojae strains on Rps1b soybeans were identified, with one represented by P7076 that carries Avr1b allele significantly changed in the coding sequence, and the second represented by P6497 that carries the intact coding sequence but deficient with Avr1b transcripts. Accumulating evidence indicate that sRNA -mediated post-transcriptional gene silencing is important for the epigenetic variation in Phytophthora. Therefore, to examine underlying mechanisms in the diploid P. sojae for the loss of Avr1b transcript accumulation in the virulent P6497, we performed sRNA deep sequencing and investigated the sRNA reads homologous to Avr1b. The results showed that there were two classes of sRNAs homologous to Avr1b (Fig. 1a, b), with one species being peaked at 25 nt (Fig. 1c) and designated as Avr1b-siRNAs. The second class was the antisense sRNA that spans the two neighboring Avr1b-siRNAs, which was named Avr1b-asRNA hereafter (Fig. 1a, b). Taken together, we infer that the silencing of Avr1b in the virulent strain P6497 may be related to these Avr1b-associated sRNA.

Fig. 1

Identification of P. sojae Avr1b sRNAs. a The abundance distribution of Avr1b sRNA reads homologous to the ORF region. The sense sRNA reads (above the x-axis) are highlighted in red, and the antisense sRNA reads (below the x-axis) are highlighted in blue. b The position distribution of Avr1b reads homologous to the ORF region. The red arrows indicate Avr1b-siRNA molecules, and the blue arrow that spans two red reads is Avr1b-asRNA. The sRNA reads and reference gene sequences are all from P. sojae strain P6497. The position and length of sRNA reads are drawn in scale. c The size distribution of Avr1b-siRNAs

Avr1b-associated sRNAs accumulate in virulent strains

To examine whether the deep sequencing identified sRNAs were associated with the regulation of Avr1b phenotype, soybean plants (rps1b) infected with P. sojae strains P6497, P7064, and P7076 were prepared and used for Northern analysis, respectively. The results showed that the Avr1b-siRNAs were specifically accumulated in the virulent strains P6497 and P7076, but not in the avirulent strain P7064 (Fig. 2a). By using the probe designed for Avr1b-asRNA, two different sizes of sRNAs were detected (Fig. 2b), with one being the Avr1b-asRNA but longer (about 40 nt) than the deep sequencing reads which ranged from 18 to 30 nt (Fig. 2b). The other one is a soybean sRNA (designated as Gma-1b-asRNA) that found in all the infection samples assayed (Fig. 2b). Avr1b-asRNA was also specifically accumulated in the virulent strains P6497 and P7076, but not in the avirulent strain P7064 (Fig. 2b). In addition, we found that Avr1b-asRNA is single-stranded, since its antisense strand was not detectable (Fig. 2c). Therefore, both Avr1b-siRNA and Avr1b-asRNA accumulation were coincided with the virulent Avr1b⁻ phenotype.

Fig. 2

Northern hybridization analyses of P. sojae Avr1b and Avr1k sRNAs during soybean infection. a The accumulation of Avr1b-siRNAs is associated with Avr1b⁻ virulence phenotype. The black triangle indicates Avr1b-siRNA. b The accumulation of Avr1b-asRNA is associated with Avr1b⁻ virulence phenotype. Except for Avr1b-asRNA (indicated by the white triangle), a constitutive host sRNA (Gma_1b-asRNA, indicated by the gray triangle) was also detected. c Northern blot failed to detect Avr1b-asRNA antisense. d Avr1k-sRNA accumulation is associated with its virulence phenotype. e Northern blot of U6 which served as a loading control. The total RNAs of 0–84 h infected tissues of G. max cultivar Williams (rps1b and rsps1k) were used. For the sample of 0 h, no P. sojae mycelia were applied. “Avr” and “avr” in the upper panel indicate avirulent or virulent genes, respectively. Each lane was loaded with 20 μg total RNA

In further, we noticed that Avr1b-asRNA was accumulated earlier than Avr1b-siRNAs during the infection (Fig. 2a, b), suggesting that Avr1b-asRNA might play as the trigger of Avr1b silencing. The detection of siRNAs homologous to Avr1b indicated that Avr1b was transcribed but silenced in the virulent strains P6497 and P7076. Thus, this species of siRNAs was likely the posttranscriptional silencing signal of Avr1b in P6497 and P7076. Compared to the soybean sRNA Gma-1b-asRNA that was instantly detected during soybean infection (or U6 in Fig. 2e), the intensity of Avr1b-siRNAs in P6497 was much stronger than that in P7076 (Fig. 2b). This probably explains why the transcripts of virulent Avr1b allele was detected easily in P7076 (Shan et al. 2004).

Expression of multiple avirulence genes in P. sojae lacks homologous sRNA

Avr1k is an avirulence gene in P. sojae strains P6497 and P7064 but virulent in strain P7076 (Song et al. 2013). We hence performed Northern hybridization to verify if sRNAs were associated with the virulent allele of Avr1k in P7076. Consistently, we found that while the homologs sRNAs were detected for the virulent Avr1k^P7076, no homologous sRNA accumulation was detected for the avirulent strains carrying Avr1k^P6497 and Avr1k^P7064, respectively (Fig. 2d). Therefore, it seems that the expression of both Avr1b and Avr1k is associated with the absence of sRNAs homologous to the avirulence gene. Inspired by these phenomena, we further examined the presence and absence of siRNA homologous to more avirulence genes. P. sojae strain P6497 carries seven cloned Avr genes. While Avr1b^P6497 is a virulent allele, all the rest six Avr^P6497 genes are avirulent (Fig. 3). The results showed that five of the six Avr genes in P6497, Avr3a/5^P6497, Avr3b^P6497, Avr3c^P6497, Avr4/6^P6497, and Avr1k^P6497, except Avr1a^P6497, had no siRNAs homologous to their open reading frames (ORFs) (Fig. 3). These results suggest that sRNAs are very likely involved in the regulation of the RXLR effector gene expressions in P. sojae.

Fig. 3

The abundance distribution of sRNA reads homologous to the functional Avr genes in P. sojae strain P6497. The sense sRNA reads (above the x-axis) are shown in red, and the antisense sRNA reads (below the x-axis) are shown in blue. The shaded regions indicate the 5′ or 3′ 500 bp flanking sequences

Bidirectional transcription mediates the production of Avr1b-asRNA

To investigate where the Avr1b-asRNA was originated from, we performed a sequence similarity search in the P. sojae genome. The result indicated that Avr1b-asRNA was derived from Avr1b or Avh1b antisense strand in the virulent strain P6497. However, in strain P7076 which lacks Avh1b (Shan et al. 2004), Avr1b-asRNA was also accumulated (Fig. 2b), indicating that the Avr1b-asRNA was derived from the antisense strand of Avr1b.

Taking together of the results above, the Avr1b region is obviously transcribed in a bidirectional manner. Sequence analysis showed that the Avr1b sense strand encodes Avr1b elicitor, while its antisense strand is predicted to encode an unknown protein with 100 amino acid residues (Fig. 4a). Hence, the dsRNA, predicted to be formed by natural sense and antisense transcripts (NATs) of Avr1b locus, is likely the precursor of Avr1b-asRNA. Avr1b-asRNA is located in the 3′ untranslated regions (UTR) of antisense Avr1b. We infer that Avr1b-asRNA is derived from the 3′-UTR part of antisense Avr1b in the NATs formed by Avr1b(+) and Avr1b(−).

Fig. 4

Bidirectional transcription of P. sojae Avr1b locus. a The model and method for detection of the bidirectional transcription of the Avr1b locus. Avr1b(+) and Avr1b(−) are organized in a tail-to-tail manner with more than 200 bp overlapping sequence. Oligo dT with an adaptor sequence (in red) was used for the first strand cDNA synthesis, and PCR was performed by using the combinations of Avr1b-asRNA sense (small blue arrow) primer, Avr1b-asRNA antisense (small black arrow) and the adaptor (small red arrow) primers. P, the promoter. b RT-PCR validation of the bidirectional transcription. G. max cultivar Williams (rps1b) infected with different strains of P. sojae were sampled at 0, 1, 2, and 3 dpi, respectively. Equal amount of total RNAs were used for each RT-PCR assays. RT(+) indicates the amplicon of RT-PCR, and RT(−) indicates the control reactions without cDNA synthesis (no reverse transcriptase applied). PsActA was amplified for 30 cycles and was used as a reference gene, while Avr1b(+) and Avr1b(−) were amplified for 35 cycles

By using strand-specific RT-PCR assay (Fig. 4a), we further verified the bidirectional transcription of Avr1b locus in P. sojae strains P6497, P7064, and P7076. We found that both the sense and antisense strands of Avr1b were transcribed (Fig. 4b). Consistently, the Avr1b-asRNA-deriving Avr1b(−) was transcribed in virulent strains P6497 and P7076 during infection, while no Avr1b(−) transcripts were detected in the avirulent strain P7064 except at the very late stage of infection (Fig. 4b). In parallel, we analyzed the expression of Avr1b(+) and there was nearly no Avr1b(+) transcripts detected in the virulent strain P6497, consistent with the results described previously (Shan et al., 2004). Both the avirulent allele Avr1b^P7064 and the virulent allele Avr1b^P7076 were detected during infection (Fig. 4b).

We conclude that the Avr1b locus is bidirectionally transcribed, leading to the development of dsRNA which may be the precursor of Avr1b-asRNA. Differential expression of Avr1b(−) is responsible for differential accumulation of Avr1b-asRNA in P. sojae (Fig. 2b and Fig. 4b). High level of Avr1b-asRNA^P6497 accumulation leads to the efficient degradation of Avr1b^P6497 transcripts and thus escaped from being detected (Fig. 4b), and Avr1b-siRNAs^P6497 are the silencing signal molecules.

Expression of large number of RXLR effector genes is associated with sRNAs

To examine whether sRNA-associated silencing of RXLR effector genes was limited to the Avr effector genes, we mapped the deep sequencing sRNAs to the entire RXLR effector genes in P. sojae. In total, 125 RXLR effector genes (about 31%) were found to be associated with the sRNA reads in strain P6497. Of these RXLR effector genes, 14 to 32 have comparable (or higher) sRNA levels than that of Avr1b sRNAs in P6497 (Fig. 5a, b).

Fig. 5

The expression of numerous RXLR effector genes was negatively associated with the accumulation of homologous sRNAs in P. sojae. Screening of the high confident sRNA-associated RXLR effector genes by sense sRNA level (a) and antisense sRNA level (b). In a and b, RXLR effector genes having comparable (or more) sRNA reads than the mean (x-axis) or max (y-axis) number of Avr1b sRNA reads homologous to the ORF region were regarded as high confident sRNA associated. Blue or green dots indicate the RXLR effector genes above the threshold. Dashed lines indicate the threshold levels for the mean or max number of reads, and the red circle indicates Avr1b. c The sRNA expression was associated with silencing of RXLR effector genes in P. sojae in mycelia and during infection. The average abundance of sense and antisense sRNA reads of the RXLR effector genes, and the expression level of the RXLR effector genes in mycelia and infection tissues (24 h post inoculation) which retrieved from Phytophthora transcriptome database were used for clustering analysis

To investigate if RXLR effector genes homologous sRNAs are associated with their silencing, we compared the expression levels of sRNA and the homologous RXLR effector genes. The results showed that the expressions of RXLR effector genes were negatively correlated with the number of sRNA reads (Fig. 5c). In particular, when the homologous sRNAs were highly expressed, nearly no transcripts of RXLR genes were detected (Fig. 5c).

Sequence variation in the regulatory regions is associated with bidirectional transcription of RXLR effector genes

Avr1b transcripts were not accumulated in virulent strain P6497, nor its antisense transcripts in the avirulent strain P7064. To examine underlying mechanism for the loss of antisense transcription of Avr1b in strain P7064, we compared the Avr1b flanking sequences in strains P6497, P7064, P7074 and P7076 that represent four major genotypes of P. sojae (Fig. 6a). The major difference within the 500-bp flanking sequences among these four strains is a 10-bp deletion in the antisense Avr1b promoter region in P7064 (deletion a, Fig. 6a). Interestingly, we also identified a 10-bp deletion in the promoter region of Avr1b in the virulent strain P7074 (deletion b, Fig. 6a), which has no (or much less) accumulation of Avr1b transcripts (Shan et al. 2004).

Fig. 6

INDEL variants identified in the regulatory regions of P. sojae genes. a Comparison of the sequence variation of regulatory regions for the Avr1b(+) and Avr1b(−). Variations in the Avr1b(+) and upstream sequences are marked in red, while that of Avr1b(−) and upstream sequences are marked in green. Avr1b-asRNA is marked in blue. R-D, the nucleotide sequence that encodes RXLR-DEER motif. b Frequency of the INDEL variants in the promoter region of P. sojae genes. The variants were called from the re-sequenced genomes of P. sojae strains P7064, P7074, and P7076, while the genome of strain P6497 was used as a reference. The 500 bp upstream sequences were regarded as the promoters for each gene. c The 9–10 bp INDEL variants are significantly enriched in the promoters of RXLR effector genes. One-sided Fisher’s exact test is used to access the significance

The two independent deletions detected in the regulatory regions suggest that these 10-bp deletions may play an important cis-regulatory role in the expression of Avr1b. This further suggests possible wide-spread presence of these 10-bp deletions in the gene regulatory regions across the genome of P. sojae. We hence performed genome-wide identification for the INDEL (insertion and deletion) variations in all the gene regulatory regions of P. sojae genes. By comparison of those INDEL variants identified in the three re-sequenced P. sojae genomes (Wang et al. 2011), we found that most of the INDELs in the promoters are 1 bp in length, and the frequency declines as the INDEL length increases (Fig. 6b). In contrast, the frequency of 9-bp and 10-bp INDEL variants were significantly increased (Fig. 6b, in red), especially in the RXLR effector genes (Fig. 6c). This suggests that the 9-bp and 10-bp INDEL variants might have undergone evolutionary selection and function as cis-elements that regulate gene expression.

Discussion

Virulence variation is critical for pathogen survival and is responsible for loss of disease resistance in a wide range of host plants against diverse pathogens. In this study, we showed that virulence variation of P. sojae is associated with accumulation of sRNAs homologous to avirulence genes. While the complete Avr1b gene sequences are present in P. sojae strains P6497 and P7064, but no Avr1b mRNA accumulation was detected during infection by P6497 (Shan et al. 2004). The Avr1b-asRNA is derived from the antisense of Avr1b, which has distinct expression patterns in strains P6497 and P7064. In the avirulent P7064, the antisense of Avr1b was not detectable in the early stage of infection. Further analysis revealed that this was likely due to a 10-bp deletion in the promoter region of Avr1b(−).

We propose a sRNA-based model (Fig. 7) to explain the results presented in this study. During the early stage of plant infection, Avr1b^P6497 transcript and its antisense transcript formed NATs, which is recognized by endogenous RNase and processed to form Avr1b-asRNA. This sRNA could lead Avr1b silencing in strain P6497 at posttranscriptional level. In the avirulent strain P7064, the lack of Avr1b(−)^P7064 leads to a failure of dsRNA formation, and therefore is unable to produce the corresponding sRNA and the silencing signals in the avirulent strain. A similar model has been proposed in Arabidopsis infected by bacterial pathogen Pseudomonas syringae carrying effector avrRpt2, in which the silencing of PPRL gene was caused by nat-siRNAATGB2, a sRNA derived from the NATs of ATGB2 and PPRL in Arabidopsis (Katiyar-Agarwal et al. 2006). In Arabidopsis, the biogenesis of siRNAATGB2 requires DCL1, RDR6, HEN1, SGS3, HYL3, and Pol IVa (Katiyar-Agarwal et al. 2006). It is quite possible that their homologs in P. sojae are also involved in the biogenesis of Avr1b-asRNA.

Fig. 7

A proposed molecular model for the regulation of Avr1b expression in the avirulent and virulent P. sojae strains. Both P. sojae strains P6497 and P7064 encodes Avr1b(+)and Avr1b(−) genes. In the virulent strain P6497 (the left panel), Avr1b(+)^P6497 and Avr1b(−)^P6497 are transcribed in a bidirectional manner and formed NATs during infection; the NATs are recognized by the putative endogenous cleavage machinery to produce Avr1b-asRNA^P6497. The produced sRNA targets on the newly transcribed Avr1b(+)^P6497 and generates abundant Avr1b-siRNA^P6497 that function as the RNA silencing signals, leading to abolished Avr1b(+)^P6497 function and escaped Rps1b recognition. In the avirulent strain P7064 (the right panel), a 10-bp deletion in Avr1b(−)^P7064 promoter region leads to significantly reduced Avr1b(−)^P7064 transcription and a failure of Avr1b(+)^P7064 silencing, allowing Avr1b(+)^P7064 be recognized by Rps1b

We have generated several lines of evidence that support our hypothesis (Fig. 7). Firstly, P. sojae encodes sRNAs that are associated with Avr1b silencing and virulence variation. By deep sequencing data analysis, we discovered two species of sRNAs: the trigger and silencing signal sRNA molecules. The Avr1b-asRNA is about 40 nt in size, different from the miRNAs and siRNAs identified in oomycetes (Vetukuri et al. 2012; Fahlgren et al. 2013; Jia et al. 2017). Although most oomycete sRNAs are 21 nt or 25/26 nt in size, larger sRNAs are also frequently observed in P. sojae and P. infestans (Vetukuri et al. 2011b; Vetukuri et al. 2012; Wang et al. 2016), and Arabidopsis (Katiyar-Agarwal et al. 2007). The siRNAs are 25–26 nt in oomycete and are recognized as the signals of gene silencing (Vetukuri et al. 2012; Qutob et al. 2013; Jia et al. 2017). Consistently, the Avr1b-siRNAs are 25 nt and associated with the virulence variation in P. sojae. In addition, we showed that a large number of sRNAs are associated with the RXLR effector gene silencing in P. sojae, suggesting that sRNA may have been associated broadly with expression regulation of RXLR effector genes. In Phytophthora, it is also shown that many sRNAs are homologous to RXLR and CRN effector genes (Vetukuri et al. 2012; Jia et al. 2017). Based on the accumulation and mapping patterns, we conclude that Avr1b-asRNA is likely the trigger sRNA and is responsible for the production of the signal siRNAs.

Secondly, bidirectional transcription and dosage compensation are simultaneously involved in the silencing of Avr1b gene. We showed that Avr1b-asRNA is derived from Avr1b(−), and thus discovered that the Avr1b locus is transcribed bidirectionally, as revealed by sequence and strand-specific RT-PCR analyses. The NATs formed by the bidirectional transcripts from Avr1b locus were the precursor of Avr1b-asRNA. NATs are commonly present and have regulatory roles in eukaryotes (Jen et al. 2005; Wang et al. 2005; Katiyar-Agarwal et al. 2006). In P. sojae, our preliminary investigation revealed that about 400 loci, including RXLR effector genes, have the potential to form NATs. Avr1b(−)^P6497 and Avr1b(−)^P7076 were abundantly expressed, being positively correlated with the supressed expression of Avr1b(+)^P6497 and Avr1b(+)^P7076. Furthermore, the expression level of Avr1b(−)^P6497 was much higher than Avr1b(−)^P7076, being consistent with the fact that Avr1b(+)^P6497 has no accumulation while Avr1b(+)^P7076 is still detectable. Consistently, Avr1b-asRNA^P6497 was much higher than Avr1b-asRNA^P7076. Therefore, it is likely that there is an expression threshold for Avr1b(−) to generate sufficient sRNA to trigger silencing of Avr1b(+).

Finally, expression of Avr1b sense and antisense transcripts are likely programmed by the INDEL variations in their promoter sequences. Our analysis of the flanking sequences revealed that the differential Avr1b(+/−) expression was associated with the INDEL variations of promoter sequences. By comparing the sequence of strains P7074 and P7064, we found that the loss of Avr1b(+)^P7074 expression is likely associated with a 10-base deletion in Avr1b(+)^P7074 promoter. By comparing the sequence of strains P7064 and P6497, the loss of Avr1b(−)^P7064 expression is likely associated with a 10-bp deletion in Avr1b(−)^P7064 promoter. A further genome-wide analysis suggests that 9 or 10 -bp INDEL variations were broadly present in the promoter regions of P. sojae genes and were significantly enriched in the RXLR effector genes, suggesting that the oomycete pathogens may have evolved functional cis-elements by promoter INDEL variations to control gene expression, especially for the RXLR effector genes. Consistently, a recent study showed that a natural 10-bp INDEL variant of the PAX7 promoter could regulate gene expression by altering the binding of transcriptional factor ZNF219 in animal cells (Xu et al. 2018). Notably, 10 bp is the distance of the double helix per turn in DNA, allowing the distance of cis-elements be easily reprogrammed by INDEL variations. Since the variation of expression of RXLR effector genes is essential for the virulence variation, these INDEL in the promoter regions may be important in shaping the plasticity of the expressions of RXLR effector genes.

Methods

Phytophthora sojae strains and culture conditions

P. sojae strains P6497 (Race 2, Avr1b⁻), P7064 (Race 7, Avr1b⁺), P7074 (Race 17, Avr1b⁻) and P7076 (Race 19, Avr1b⁻) representing four major genotypes and nearly all the genetic variations (Forster et al. 1994) were used in this study. P. sojae cultures were routinely maintained on 5% (v/v) carrot juice agar (CA) medium supplemented with 0.002% (w/v) β-sitosterol and 0.01% (w/v) CaCO₃, and incubated at 25 °C in the dark.

Soybean inoculation procedure

Soybean (Glycine max) cultivars Williams (rps1b) and Williams L77–1863 (Rps1b) were used for P. sojae infection assay. The 8-day-old seedlings were wound inoculated at hypocotyls with a vertical slit of about 1 cm in length, and a small piece of P. sojae mycelia were then mounted at the wound. The hypocotyls were cut from the upper 0.5 cm to the lower 0.5 cm (total 2 cm including the wound) 0 to 3 days post inoculation (dpi). Williams infected tissues were used for Northern blot and RT-PCR analyses, while Williams L77–1863 carrying Rps1b resistance gene was used as a control to monitor the phenotype.

Phytophthora sojae and soybean RNA preparation

RNA was isolated by using TRIzol (Invitrogen, USA). In brief, the infected hypocotyl tissues were ground into fine powder with liquid nitrogen, treated with TRIzol at room temperature, followed by extraction with chloroform and phenol-chloroform. RNA precipitation was performed by isopropyl alcohol overnight at − 20 °C to fully recover sRNAs. Finally, the sRNA pellets were washed with ethanol and were dissolved with DEPC-treated ddH₂O.

Northern assay for sRNA and siRNA detections

Northern analysis was performed as described previously (Wang et al. 2016). In brief, 20 μg of total RNAs were denatured and subjected to a 15% urea-denaturing PAGE. RNAs were transferred onto nylon membranes and were fixed by a UV-crosslinker (HL-2000 HybriLinker, UVP) followed by baking at 80 °C for 2 h. Pre-hybridization and hybridization were carried out at 42 °C for 8 h, respectively, in a Hybridizer (HL-2000 HybriLinker, UVP). The membranes were washed twice by 2 × SSC and 0.5% SDS at 37 °C for 15 min. The images were exposed and detected by the FLA-7000 image system (Fujifilm, Japan).

For Avr1b-asRNA detection, DNA oligonucleotides (Avr1b-asRNA-probe: 5′-TCA GTC ACG CTG AAG GTC CGC TCC C-3′), complementary to the Avr1b-asRNA, was labeled with [γ-³²P]-ATP using T4 Polynucleotide Kinase (TaKaRa, China). For Avr1b-siRNAs detection, the ORF sequences were amplified from the genomic DNA and labeled with [α-³²P]-dCTP using Random Primer DNA Labeling Kit (TaKaRa, China).

RT-PCR analysis

RT-PCR assay was performed to verify gene transcription. To confirm the expression of Avr1b and its antisense transcript, the first strand cDNA was synthesized using 1 μg total RNA with oligo dT carrying an adaptor (OldTail: 5′-GAG CAG TCG AGA CTC GAT CGC ATG CAT TTT TTT TTT TTT TTT TTV N-3′) by using Reverse Transcriptase M-MLV (TaKaRa, China). Two gene-specific primers (GSP) designed according to the Avr1b and its antisense sequences (GSP+: 5′-GGG AGC GGA CCT TCA GCG TGA CTG A-3′; GSP-: 5′-TCA GTC ACG CTG AAG GTC CGC TCC C-3′) and the adapter primer (Tail: 5′-GCA GTC GAG ACT CGA TCG CAT-3′) were used as the forward and reverse primers, respectively.

Bioinformatics analysis

For the identification of Avr1b-asRNA and Avr1b-siRNA, the sRNA reads (SRR4450568) (Wang et al. 2016) were mapped to the open reading frame (ORF) of Avr1b^P6497 by using bowtie (Langmead et al. 2009) with no mismatch allowed as described previously. Mapping of the reads to the entire RXLR effector genes was done as that of Avr1b. The expression data of RXLR effector genes were retrieved from Phytophthora transcriptome database (Ye et al. 2011).

For the analysis of Avr1b regulatory sequences, Avr1b^P6497 and its flanking sequences were retrieved from P. sojae genome (P6497, version 3) (Tyler et al. 2006) and were used for identification of the contigs in P. sojae strains P7064, P7074, and P7076 (Wang et al. 2011). The flanking sequences of P7064 and P7074 were derived from single contigs (contig29872 and contig03312, respectively) and P7076 from two separate contigs (contig01119 and contig02466) in the 454 assemblies. The retrieved sequences were confirmed by sequencing at GenScript (Nanjing, China). Multiple sequence alignment was performed with Clustal X2 (Larkin et al. 2007) and visualized via DOG (Ren et al. 2009) The promoter sequences were predicted by EP3 (Abeel et al. 2008).

The identification of INDEL variants was performed as described previously (Wang et al. 2017). Briefly, the contigs of the re-sequenced P. sojae genomes of strains P7064, P7074, and P7076 (Wang et al. 2011) were re-ordered by ABACAS (Assefa et al. 2009) with the reference genome of strain P6497 (Tyler et al. 2006), and MUMmer (Delcher et al. 1999) was used to call the INDEL variants with its default parameters. According to previous study, the 5′ UTR regions are very short (Win et al. 2006), thus the 500-bp upstream sequences were regarded as the promoters for each gene.