Phytophthora sojae strain P6497 was maintained on 5% CA medium at 25°C in the dark. To prepare P. sojae-infected tissues, soybean (Glycine max cv Williams) seedlings of 7–10 days were wound-inoculated with P. sojae mycelia at hypocotyls as described (Tyler et al., 1995; Shan et al., 2004a). The infected hypocotyls of 48 hpi (hours post inoculation) were collected for expression pattern analysis.

TRIzol (Invitrogen, USA) reagent was used for total RNA isolations following the standard protocol except that RNA precipitation was carried out at –20°C overnight.

Equal amounts of RNA from P. sojae tissues representing five distinct developmental stages, including oospores, mycelia, sporangia, cysts, and germinated cysts were mixed and sequenced at Beijing Genomics Institute (18–30 nt, accession number: SRR4450568). Sequencing of the 18–50 nt sRNA, the 28–50 nt sRNA, and the degradome of P. sojae mycelia were performed at LC Sciences (accession number: SRR4450567, SRR4450570, SRR4450569, respectively).

A total of 1,360,211 clean P. sojae sRNA reads were obtained by deep sequencing. All reads aligned to the P. sojae genome (Tyler et al., 2006) were used for further analysis. Ribosomal RNAs were retrieved from GenBank and gene models generated by the P. sojae genome project were also used to filter the short tags. The tsRNAs were identified as shown in Supplemental Figure S1. In brief, eukaryotic tRNA sequences downloaded from Genomic tRNA Database (Chan and Lowe, 2009) were used to extract tRNA-derived sRNA reads. The 75 nt upstream and downstream flanking genomic sequences were extracted for clover-leaf structure screening by RNAshapes (Steffen et al., 2006). We chose the highest MFE (minimum free energy)/length ratio fragment of the corresponding genomic sequences and assembled them to unique sequences by CAP3 (Huang and Madan, 1999) to generate a reference sequences set of tRNA-derived sRNA tags, and the precise tRNA coding regions were defined by feeding the non-redundant reference set to tRNAscan-SE (Lowe and Eddy, 1997). The sRNA sequences were mapped to their references using Bowtie (Langmead et al., 2009) with no mismatch allowed. The generated landscapes were used to determine the sRNA candidates by identification of the peaks in the mapping landscapes via the Base Deep Index (BDI = floor[9×(x-x_min)/(x_max–x_min)]) we defined, which reflects the relative coverage of each base along the tRNA sequence. All tsRNA core sequences identified were listed in Supplemental Table S1.

To search for potential target genes of tsRNAs in P. sojae, the entire gene models with their 200 bp flanking sequences were searched against tsRNA candidates using NCBI BLAST 2.2.18 (Altschul et al., 1997). Since no information was available on the presence of sequence-specific targeting by tsRNAs in Phytophthora, we used relatively less stringent criteria and searched potential target transcripts that have complementary match in local perfect manner with tsRNAs. Only antisense hits (E-value = 10, leading at least 12 bp) were considered as candidate targets. The details of tsRNA-target pairs and the potential function of target genes is shown in Supplemental Tables S2 and S3, respectively. The relationship between tsRNAs and their targets were visualized by Circos (Krzywinski et al., 2009). Digital transcriptome data were retrieved from the Phytophthora Transcriptome Database (Ye et al., 2011) and were used to evaluate the expression patterns of P. sojae tsRNA targets. Target genes without expression data in the five developmental stages were excluded from the analysis, and the resulting matrix was normalized by both row and column to generate the heatmap using Pearson’s correlation coefficient for pairwise distances. The selected 10 targets for qPCR verifications are listed in Supplemental Table S4. The degradome sequencing reads were mapped to the reference sequence via Bowtie (Langmead et al., 2009).

The re-sequenced P. sojae tsRNAs were analyzed using in-house tsRNA annotation pipeline tsRFinder¹ with its default parameters. All the reference tsRNA sequences were list in Supplemental Table S5.

For manual cloning of tsRNAs, sRNAs ranging from 15 to 50 nt in size were isolated from the total RNA extracted from mixed P. sojae oospores, mycelia, sporangia, cysts, and germinated cysts. After dephosphorylation with Antarctic Phosphatase (NEB, USA), the sRNAs were ligated with 3′ AppDNA adaptor A1 by T4 RNA ligase (NEB, USA) at 14°C overnight, and separated by 15% Urea-PAGE. T4 Polynucleotide Kinase (NEB, USA) was used to add additional 5′-phosphates to the recycled products before subsequent ligation with 5′ RNA adaptor A2 by T4 RNA ligase (NEB, USA). The ligated sRNAs were reverse transcribed with M-MLV (Promega, USA) using 3′ PCR primer P1. The reversely transcribed sRNAs were amplified with primers P1 and P2 and the products were cloned into pGXT vector (Chen et al., 2009) and sequenced at Beijing Genomics Institute. The sRNAs homologous to the tRNAs were regarded as tsRNAs. The adaptors and primers used for sRNA cloning are listed in Supplemental Table S6.

To perform Northern analysis, 20 μg of total RNAs were denatured at 80°C for 10 min and then subjected to gel separation in a 15% urea-denaturing PAGE. sRNAs were electro-blotted onto Amersham Hybond-N⁺ membrane (GE Healthcare, UK) with a semi-dry transfer cell (BioRad, USA) at 3.3 mA/cm² for 20 min. RNA was fixed by a UV-crosslinker (HL-2000 HybriLinker, UVP) at 1200 × 100 μJ/cm² energy for 2 min, followed by baking at 80°C for 2 h. After pre-hybridization with 10 mL hybridization buffer consisting of 1.5 mg Herring Sperm DNA, 7% SDS and 0.2 M Na₂HPO₄ for 8 h, DNA oligonucleotides complementary to the tsRNAs (Supplemental Tables S1 and S5) were labeled with [γ-³²P]-ATP using T4 Polynucleotide Kinase (TaKaRa, China) and added to the hybridization solution. Hybridization was carried out at 42°C for 8 h in a Hybridizer (HL-2000 HybriLinker, UVP). After discarding the hybridization solution, the filter blots were washed twice, each with 30 mL 2X SSC and 0.5% SDS (preheated at 37°C) for 15 min at 37°C. The filter blots were exposed to phosphoimager screens and the signals were detected by scanning with the FLA-7000 image system (Fujifilm, Japan).

For real time RT-PCR, 500 ng of total RNA was treated with DNase I (TaKaRa, China) to remove genomic DNA and the cDNA synthesis was performed using PrimeScript^TM RT reagent kit according to the manufacturer’s instructions. PCR reactions were performed using 5 μL of a 1:20 dilution of the first-strand cDNA with SYBR Premix Ex Taq^TM II according to the manufacturer’s instructions (TaKaRa, China). Expression levels of the tested genes were quantified using an iQ5 Real-Time Cycler (Bio-Rad, USA). The P. sojae gene PsActA (gene id 108986) was chosen as a reference gene (Zhao et al., 2011). All samples were standardized to the levels in mycelia and the fold changes of expression were calculated using the 2^ΔΔCt method. The data were calculated or analyzed by customized R scripts. All the primers used for quantitative RT-PCR analysis are listed in Supplemental Table S6.

To perform 5′ RNA ligase -mediated rapid amplification of cDNA ends (5′ RLM-RACE), the poly (A)⁺ RNA was isolated by using Magnetic mRNA Isolation Kit (NEB, UK) from 200 μg total RNA. The cleaved products are uncapped and carry free phosphate, thus enabling direct ligation with an RNA adaptor RA44 by using T4 RNA Ligase (Ambion, USA). The ligation products were extracted by phenol/chloroform and precipitated with glycogen before synthesis of the first-strand of cDNA by using SuperScript II Reverse Transcriptase (Invitrogen, USA). A nested PCR with ExTaq Hot Start Version (TaKaRa, China) was used to detect the degraded products with primers RA44OP/IP and GSP1/GSP2. The amplicons were further confirmed by sequencing. All adaptors and primers used for 5′ RLM-RACE analysis can be found in Supplemental Table S6.

For TS4WT, TS4MU, TS7WT, and TS7MU, the GUS gene was amplified using primer pairs GUSF/GUSR (carrying restriction enzyme sites and tsRNA-binding sites, see Supplemental Table S6) from pCXGUS-P (Chen et al., 2009), double digested with Xmal I and Nsi I (Promega, USA), and inserted into plasmid pUbIN. The GUS-BS fragments were then released by Xma l and Nsi I, blunted with Pfu DNA polymerase (TaKaRa, China), and inserted into the Sma I site of pHAM34H. For TS7WT2 and TS7MU2, the GUS gene was amplified and inserted directly into the Sma I site of pHAM34H. These GUS sensor expression constructs were transiently expressed in P. sojae protoplasts with polyethylene glycol as described (Dou et al., 2008; Fang and Tyler, 2016). Up to 30 μg test plasmid and 10 μg selection marker plasmid were used. GUS activities were determined at 12 h post incubation using histochemical GUS staining (Dou et al., 2008).

To identify Phytophthora sojae tsRNAs, we surveyed sRNAs from five distinct developmental stages of the P. sojae lifecycle, including sexual oospores, asexual mycelia, sporangia, cysts, and germinated cysts. High throughput sequencing led to the identification (see Supplemental Figure S1) of a total of 41 candidate tsRNAs derived from 26 tRNAs (Figure 1A and Supplemental Table S1). To confirm the existence of tsRNAs, we performed small-scale manual cloning of tsRNAs in the 15–50 nt RNA fractions (Figures 1B,C; Supplemental Table S1). We chose one tsRNA identified by deep-sequencing and one identified by manual cloning to further confirm existence of P. sojae tsRNAs. The results showed that both the cloned and sequenced tsRNAs were detected by Northern blots (Figure 1D). These three independent assays indicate that the tRNA-derived fragments we detected are real endogenous sRNAs in P. sojae. Northern analyses also indicate that the two tsRNAs selected here are conserved in another oomycete pathogen, P. parastica (Figure 1D). Aligning the tsRNA sequences to their source tRNA sequences indicated that these tsRNAs were derived from the 5′ arm, the 3′ arm or both arms of the tRNA (Figure 1E; Supplemental Table S1). We named the tsRNAs by their source tRNAs, for instance, tsRNA-IleAAT is derived from tRNA-IleAAT.

To explore the potential biological roles of tsRNAs, we examined the accumulation of tsRNAs at six different developmental stages of the P. sojae lifecycle. The result showed that these tsRNAs had differential levels of accumulation at the different stages (Figure 2). For example, tsRNA-IleAAT displayed relatively high levels of accumulation in oospores, mycelia, and sporangia, nearly undetectable levels in cysts and in germinated cysts, low levels upon plant infection (Figure 2A, left panel) with gradually increasing levels during plant infection (Figure 2A, right panel). Strikingly, such expression patterns were nearly identical for most of P. sojae tsRNAs during these development stages (Figure 2B, Supplemental Figure S2 and Table S1). The expression levels of tsRNAs relative to cognate tRNAs were highly variable, suggesting that the relative amount of tsRNAs is independent of tRNA abundance (Supplemental Figure S1 and Table S1).

To identify if there are potential homologous target genes of tsRNAs, we searched all available P. sojae gene models with their 200 bp flanking sequences for complementary matches to the candidate tsRNAs. A total of 456 complementary sites (Supplemental Table S2), representing 371 putative tsRNA target genes with various functions (Supplemental Table S3), were identified. Notably, we found that one tsRNA could target to multiple genes, and one gene could bind with multiple tsRNAs (Figure 3). Among the tsRNA candidates examined in the analysis, 32 had the potential to bind to at least one transcript (Figure 3).

To investigate whether there is a negative correlation between tsRNA accumulation and target gene expression, we searched the Phytophthora Transcriptome Database (Ye et al., 2011) and obtained digital gene expression profiles of the predicted tsRNA target genes at five developmental stages of P. sojae: mycelia, sporangia, cysts, germinated cysts and infection (24-hpi soybean leaves). Since most tsRNAs have a similar and stage-specific developmental expression pattern (Figure 2B), a similar but inverse expression pattern of the putative targets was expected. Cluster analysis based on gene expression levels showed that, the biggest cluster of target genes, consisting of 70% (126/181) of the genes analyzed, exhibited abundant transcripts in cysts in which tsRNAs were nearly undetectable, and extremely low transcript levels in mycelia and sporangia where tsRNAs were abundant (Figure 4A). Thus, for the majority of the putative target genes, their expression levels correlated negatively with the abundance of tsRNAs. To confirm the profile and negative correlation, we analyzed the expression of 10 putative tsRNA target genes (Supplemental Table S4) using quantitative RT-PCR, and nine of them showed highest expression level in cysts (Figure 4B), consistent with the digital expression profile (Figure 4A). For example, tsRNA-GluCTC-2 target P. sojae gene transcript TS7, while tsRNA-GluCTC-2 is absent in cysts, its target TS7 is highly up-regulated in this stage (Supplemental Figure S3).

To investigate whether the inverse correlation between tsRNA accumulation and target transcription levels was due to tsRNA mediated mRNA degradation in P. sojae, we employed Parallel Analysis of RNA Ends (PARE) technology and sequenced RNA degradome in P. sojae mycelia, in which tsRNAs are highly accumulated. We obtained a total of 26,039 reads perfectly matched to 160 predicted tsRNA targets (43% of the total predicted), and 95 of them belong to the biggest cluster (75% of the 126 genes) which were negatively correlated with accumulation of homologous tsRNAs. Most of these 95 genes had fewer matched reads, which were likely a result of tsRNA-mediated degradation. Examining the distances between the tsRNA binding sites and the degradome reads showed that, the degradation events occurred mainly around the binding sites, but extended toward the 3′ direction (Figure 5A). Of the 10 genes selected for quantitative RT-PCR confirmation, eight of them had few reads mapped. For example, both TS4 and TS7 had a relative many more reads around the binding sites, and many reads were found at the 3′ downstream (Figure 5B), indicating that their down-regulation in mycelia is likely due to tsRNA mediated transcript degradation. To verify the degradome sequencing results, we performed 5′ RLM RACE to capture the degraded RNA products in the mycelia of P. sojae. Consistent with the degraded reads patterns, we observed that the degradation events occurred around the binding sites with a little bias to the downstream 3′ region (Supplemental Figure S4).

The degradome results suggested that tsRNAs mediate RNA degradation in a homology-dependent manner. Since tsRNA is highly expressed in mycelia and is generated from multiple copies of tRNA, it is difficult and unnecessary to knock-down and over-express tsRNA to determine the homology-dependent RNA degradation. To further examine this possibility, we generated sensor constructs containing the β-glucuronidase (GUS) reporter gene fused with either the wild-type (WT) or mutant (MU) sequences of the tsRNA-SerGCT2 and tsRNA-GluCTC2 binding sequences of target genes TS4 and TS7, respectively (Figure 6A). In the mutant sequences we introduced multiple nucleotide mismatches to disrupt base-pairing with the respective tsRNAs. We transiently transformed P. sojae protoplasts with the sensor constructs and determined GUS expression using histochemical staining. Almost no blue staining or GUS activity was observed in the independent P. sojae cell lines transformed with the TS4WT construct containing the wild-type tsRNA-binding sequence (Figure 6B). In contrast, P. sojae cell lines transformed with the mutant construct TS4MU showed clear blue staining indicating good levels of GUS expression.

We assayed two sets of GUS sensor constructs for testing tsRNA-GluCTC2-mediated gene silencing, TS7WT, TS7MU and TS7WT2 and TS7MU2. We first prepared and tested the TS7WT and TS7MU constructs containing a short version of the target sequences. The result showed that the WT construct gave very limited decrease in GUS sensor expression in comparison to the mutant construct (Figure 6B, left panel). We then considered that sequences flanking tsRNA-binding site could impact the efficiency of tsRNA-mediated gene silencing, and therefore prepared and tested the second set of constructs, TS7WT2 and TS7MU2, containing the tsRNA-binding sites plus 10-nt 5′ and 3′ flanking sequences (Figure 6A). This modification dramatically increased the sensitivity of the GUS sensors, and P. sojae cell lines transformed with the TS7WT2 construct showed strongly reduced GUS activity in comparison to those transformed with the TS7MU2 construct (Figure 6B, right panel). To determine if the repressed GUS activity in the wild-type fusion transgenes is a result of mRNA degradation, we performed quantitative RT-PCR analysis of the transformed P. sojae cells. GUS sensor mRNA levels in transformed P. sojae cells were much (∼15-fold, P-value <0.01) lower with the WT construct than the mutant construct (Figure 6C), indicating that GUS sensor repression is associated with RNA degradation. Taken together, results with the GUS sensor constructs indicate that P. sojae tsRNAs are capable of triggering RNA degradation-based gene silencing, and suggest that this tsRNA-induced gene silencing requires high-level sequence complementarity between tsRNA and target RNA.

Owing to the short-range size selection (18–30 nt) in the original sRNA sequencing (Figure 1A), the tsRNA dataset used in this study contained mostly the core sequences of tsRNAs but not the full-length size (>30 nt), as expected from the northern blot hybridization results (Figure 2; Supplemental Figure S2). To obtain an improved high-quality reference tsRNA sequences, we re-sequenced two sRNA libraries from tsRNA-enriched P. sojae mycelia, PsMa and PsMb, with selected sizes ranging from 18–50 nt and 28–50 nt (Figure 7A), respectively. In these sRNA datasets, the size of tsRNA reads ranged from 29–36 nt and peaked at 33 nt (Figure 7B). Interestingly, the 5′ terminal nucleotides of the tsRNA reads appear to be biased toward G (Figure 7B), which is different from miRNA and siRNA which tend to have U at the 5′ end. It is also notable that up to 12% of the total sRNA reads in the range of 28–45 nt were derived from tRNA (Figure 7C, right panel).

By identification of the tsRNAs in these two sRNA libraries, an improved high-quality tsRNA repository was generated (Supplemental Table S5). These tsRNAs are generally highly expressed in the mycelia of P. sojae. For example, tsR1 (Figure 8A, highlighted in blue), derived from the 5′ arm of tRNA-GlyCCC1, is the most abundant tsRNA in PsMb, accounting for 2.05% of the total reads (Supplemental Table S5). The 3′ arm of tRNA-GlyCCC1, tsR18, was also highly expressed (Figure 8A, highlighted in red) in the repository (Supplemental Table S5). To verify if the tsRNAs identified from re-sequencing are full length mature tsRNAs, we performed Northern blot analysis with tsR1. The result is consistent with the size detected by Northern blotting being identical with that detected by re-sequencing (Figure 8B).