Botrytis cinerea small RNAs are associated with tomato AGO1 and silence tomato defense-related target genes supporting cross-kingdom RNAi

Cross-kingdom or cross-species RNA interference (RNAi) is broadly present in many interacting systems between microbes/parasites and their plant and animal hosts. A recent study by Qin et al. (2022) performed correlation analysis using global sRNA- and mRNA-deep sequencing data of cultured B. cinerea and B. cinerea-infected tomato leaves and claimed that cross-kingdom RNAi may not occur in B. cinerea–tomato interaction (Qin et al., 2022). Here, we use experimental evidence and additional bioinformatics analysis of the datasets produced by Qin et al. (2022) to identify the key reasons why a discrepancy between the conclusion of Qin et al. 2022 and previously published findings occurred. We also provided additional experimental evidence to support the presence of cross-kingdom RNAi between tomato and B. cinerea. We believe it is important to clarify the basic concept and mechanism of cross-kingdom/cross-species sRNA trafficking and illustrate proper bioinformatics analyses in this regard for all the scientists and researchers in this field.


INTRODUCTION
Cross-kingdom or cross-species RNA interference (RNAi) is a phenomenon in which small RNAs (sRNAs) are transported between interacting organisms to silence the expression of target genes in the adversary or interacting partner (Weiberg et al., 2013, Cai et al., 2018, Weiberg et al., 2014. Our lab first discovered this communication mechanism between plants and the fungal pathogen Botrytis cinerea. Some sRNAs of B. cinerea are transported into host plants where they hijack host RNAi machinery for cross-kingdom gene silencing (Weiberg et al., 2013). We found that the transferred fungal sRNAs are loaded into the host Arabidopsis Argonaute 1 (AGO1) protein to silence host target genes involved in immune responses (Weiberg et al., 2013). This phenomenon has since been observed later in many different interacting organisms. For example, other plant fungal pathogens (such as Fusarium oxysporum and Verticillium dahliae), oomycete pathogens (such as Hyaloperonospora arabidopsidis), the parasitic plant Cuscuta campestris, and even the symbiotic ectomycorrhizal fungus Pisolithus microcarpus and the bacterial symbiont Rhizobium (Bradyrhizobium japonicum) can send a set of sRNAs into the host plant to silence the expression of host target genes (Dunker et al., 2020, Shahid et al., 2018, Ren et al., 2019, Wong-Bajracharya et al., 2022, Ji et al., 2021, Zhang et al., 2022. Furthermore, the mode of action of transferred microbial sRNAs seems to be conserved. The sRNAs from F. oxysporum (Ji et al., 2021), V. dahliae (Wang et al., 2016), H. arabidopsidis (Dunker et al., 2020)

and the
Rhizobium (Ren et al., 2019) also utilize host AGO proteins for silencing host target genes.
Excitingly, cross-kingdom or cross-species RNAi has also been observed in animal-pathogen or parasite interaction systems. For example, the fungal pathogen of mosquito, Beauveria bassiana transfers a microRNA-like RNA (milRNA) to the host cells. This milRNA is loaded into mosquito AGO1 to silence the host immunity gene Toll receptor ligand Spätzle 4 (Cui et al., 2019). 4 Some parasites, such as the gastrointestinal nematode Heligmosomoides polygyrus, send a range of its sRNAs into mammalian gut cells to target and silence immunity and inflammation-related genes (Buck et al., 2014, Chow et al., 2019. The causal bacterial pathogen for human Legionnaires' disease, Legionella pneumophila also sends bacterial sRNAs into human host cells to regulate the expression of innate immune response genes (Sahr et al., 2022). Thus, crosskingdom/cross-species RNAi is widely observed and accepted by scientists studying many interaction systems between microbes/pests and their plant or animal hosts.
Recently, Qin et al. (2022) performed sRNA-and mRNA-deep sequencing analysis of cultured B. cinerea and B. cinerea-infected tomato leaves and conducted inoculation assays with mutated B. cinerea strains that deleted a transposon or two Dicer-like genes (bcdcl1dcl2). They claimed that cross-kingdom RNAi does not occur in B. cinerea-tomato interaction (Qin et al., 2022). Here, we re-analyzed RNA sequencing data released by this study and aimed to recapitulate infection assay results to identify the key reasons for the discrepancy with previously published findings. We also provided additional experimental evidence that B. cinerea sRNAs are associated with tomato AGO1 to silence tomato target genes, which supports the presence of cross-kingdom RNAi between tomato and B. cinerea. We believe it is important to clarify the basic concept and mechanism of cross-kingdom/cross-species RNAi and to illustrate the proper bioinformatics analyses of sRNA and mRNA deep sequencing libraries. Among them, 7,042 B. cinerea sRNAs were predicted to target 3,185 distinct tomato mRNAs.

Discrepancies in the bioinformatics analysis and interpretation of sRNA and
Conversely, 114,011 tomato sRNAs were found that were predicted to target 11,434 B. cinerea mRNAs. Based on the sRNA-seq and mRNA-seq data, the authors observed a weak negative correlation between the entire set of sRNAs and their potential targets that was not significant according to a permutation analysis. Thus, they concluded that there was no evidence for crosskingdom RNAi in the B. cinerea-tomato interaction. This conclusion is based on the assumptions that all sequenced B. cinerea sRNAs can contribute to cross-kingdom RNAi, that an important part of the plant gene expression changes observable by mRNA-seq are caused by fungal sRNAs, and that predictions of sRNA-mRNA interactions have low levels of false positives. We believe that these assumptions are not accurate for the following reasons: First, not all sequenced sRNAs of B. cinerea are genuine extracellular sRNAs, not all extracellular sRNAs would translocate into the host cells, and not all sRNAs that enter host cells will form a functional silencing complex, e.g., with the host AGOs. Even within the same organism, not all sRNAs can silence their bona fide targets at any time. Many sRNAs, such as most mammalian miRNAs, can have an average of 300-400 predicted targets; some even have over 1000 predicted target genes. However, many studies have shown that only a small subset of these target genes' transcripts are reduced in abundance at any given time under specific developmental 6 stages or in response to various environmental cues. In contrast, most predicted target genes are not suppressed (Bartel, 2009, Wilczynska & Bushell, 2015. This phenomenon cannot lead to a conclusion that there is no evidence for mammalian miRNA to induce the silencing of target genes.
Thus, conducting global sRNA-mRNA correlation analysis for cross-kingdom/cross-species RNAi would not yield meaningful results.
Second, the quality and reproducibility of Qin et al.'s sRNA-and mRNA-seq datasets are questionable. The sRNA data includes large amounts of reads mapping to rRNA and tRNA (up to over 90%) for all their analysis (see below for details). This could be an indication of low-quality and highly degraded RNA, and these sequences will mask any real silencing patterns. Restricting the analysis to sRNAs of 20-24nt, as per Qin et al.'s methods, will help increase the signal-tonoise ratio, as would excluding some reads mapping to rRNA/tRNA. Nevertheless, these steps will not effectively remove the problem, since all the RNA in these samples would have suffered similar levels of degradation, and even 20-24nt reads would include a large fraction of irrelevant degradation products.
Third, Qin et al. predicted tomato targets using the psRNATarget tool with the default settings of V2 2017, without any further requirements. These settings allow bulges and gaps in the target site base pairing and up to two mismatches in the seed region. An Expectation score of 5 was used as a cutoff. The stringency of these settings for plant target prediction is rather low because gene silencing rarely occurs with a bulge or gap in the sRNA-target pairing regions and mismatch at the 10 th -11 th nt positions will block the function of sRNA (Schwab et al., 2005), which was the principal basis for designing the short tandem target mimic (STTM) to block the function of sRNAs (Todesco et al., 2010). To achieve more precise target prediction, in our previous paper (Weiberg et al., 2013), we used the TAPIR1.1 tool (Bonnet et al., 2010, Schwab et 7 al., 2005 to predict tomato targets with more stringent requirements: a) No gap or bulge within the alignment between the sRNA and the target was allowed; b) The 10 th and 11 th nts of the sRNA must perfectly match its target; c) At most one mismatch or two wobbles were allowed from position 2 to 12; d) A maximum of two continuous mismatches was allowed; e) A score of 4.5 was used as a cutoff.
Thus, the extremely high number of tomato target genes predicted by Qin et al. (2022) was most likely overestimated. Furthermore, many of their predicted target genes in tomato may not be true Bc-sRNA targets. Importantly, even true sRNA/mRNA targets will not be functional unless at the very least, they are present in the same cell and within a silencing complex.
Most importantly, the essential step after the deep-sequencing and bioinformatics analyses is to perform experimental validation analysis, preferably using multiple approaches, to confirm all the conclusions drawn from any omics analysis. All the studies from our labs have followed this principle, and we always use molecular, biochemical, cell biological and/or genetics approaches to confirm any results we obtained from omics data.

Some B. cinerea sRNAs associate with tomato AGO1 supporting the presence of crosskingdom RNAi between B. cinerea and tomato.
We previously observed that a subset of B. cinerea sRNAs move into plant cells and utilize Arabidopsis AGO1 proteins to silence host target genes (Weiberg et al., 2013). Following Arabidopsis AGO1 co-immunoprecipitation and RNA isolation, the B. cinerea Bc-siR3.1, Bc-siR3.2 and Bc-siR5 were detected by a stem-loop RT-PCR method (Weiberg et al., 2013). The same molecular mechanism was also found in plant pathogens F. oxysporum, V. dahliae, H. arabidopsidis, animal fungal pathogen Beauveria bassiana and even in beneficial bacterial 8 symbiont Rhizobium during cross-kingdom RNAi with their hosts (Ji et al., 2021, Zhang et al., 2022, Dunker et al., 2020, Ren et al., 2019. The sRNAs from these microbes are loaded into host AGO proteins to silence host target genes. In order to experimentally validate whether cross-kingdom RNAi is present between B. cinerea-tomato during infection, the most straightforward way is to examine whether any B. cinerea sRNAs are associated with tomato AGO proteins after moving into tomato cells by performing tomato AGO1-immunoprecipitation (AGO1-IP) analysis. We utilized the tomato AGO1 antibody generated by the Katherine Borkovich lab, which was used to successfully pulldown tomato sRNAs (Ji et al., 2021), and immunoprecipitated the sRNAs associated with the tomato AGO1 after B. cinerea infection.
We indeed detected some incorporated Bc-sRNAs. The three Bc-sRNAs described in our previous study (Weiberg et al., 2013), Bc-siR3.1, Bc-siR 3.2, and Bc-siR5, were detected to be associated with tomato AGO1 complex at 12, 16 and 24 hpi ( Figure 1a). We also checked Bc-sRNAs predicted to target tomato genes by Qin et al. (2022) and mapped them to annotated regions in the B. cinerea B05.10 genome. We found that two Bc-sRNAs (#3 and #5) were also associated with tomato AGO1 (Figure 1b). To determine whether these AGO1-associated Bc-sRNAs can trigger host target gene silencing, we examined the target gene expression after B. cinerea infection.
Silencing of a Bc-sRNA#5 target, the receptor-like kinases (SIRLK) was also detected by Qin et 9 al. in figure 4 of Qin et al. (2022). Most of these Bc-sRNA tomato target genes are associated with plant immune responses. For instance, ATG2 negatively affects powdery mildew resistance and mildew-induced cell death in Arabidopsis (Wang et al., 2011). Knockdown of SlMPKKK4 using virus-induced gene silencing (VIGS) resulted in enhanced disease susceptibility in response to B.
cinerea (Weiberg et al., 2013). The SlACIF1 is involved in Cladosporium fulvum infection in tomato (van den Burg et al., 2008). Some PPR-derived siRNAs have been predicted to target genes in Verticillium dahliae, indicating their potential contribution to defense against fungal pathogen infection (Hudzik et al., 2020). Many of the RLKs have been found to detect microbe-and hostderived molecular patterns, serving as the first layer of defense responses (Tang et al., 2017). The transcript accumulation of their predicted tomato target genes was reduced during the infection, supporting that Bc-sRNA-mediated silencing of target genes occurred in tomatoes. Thus, we demonstrated that cross-kingdom RNAi is evident between B. cinerea and tomato. However, many Bc-sRNAs identified by Qin et al. (2022) were not associated with tomato AGO1 (#4, #6, #8, #9 in Figure 1b), which could explain why they did not silence the expression of corresponding predicted targets in tomato ( Figure 1c).

The sequencing data showed poor quality and a lack of reproducibility
To fully understand the results and conclusions from Qin et al. (2022), we decided to reanalyze their deep-sequencing data regarding quality and reproducibility. In order to avoid the appearance of potential bias, we invited another two independent hard-core bioinformatics groups the first replicate of cultured B. cinerea (B16_1) shows a relatively uniform sRNA size distribution across 18-38 nt, with no obvious peak around 20-24nt. This is a typical pattern for degraded RNA, which is corroborated by the high levels of rRNA and tRNA-mapping reads. The first replicate of the ku70 sample differed from the other two ku70 replicates. In the first ku70 replicate, a significant proportion of 22 nt sRNAs originated from TEs, while the other two ku70 replicates showed a much lower proportion of sRNAs at this length. Additionally, the pattern of the ku70 replicates was noticeably distinct from that of the cultured B05.10 B. cinerea samples, further demonstrating the big difference between the B05.10 and ku70 mutant strains. The first two replicates of infected tomato samples (I12As and I12Bs), show a main peak at 33 nt, mainly consisting of tRNA fragments., while the third replicate (I12Cs) shows a substantially different profile ( Figure 2b). If we assume that TE-derived sRNAs are likely to be responsible for cross-kingdom RNAi, the proportion of such reads varies drastically between replicates. We found that the vast majority of Second, we also performed PCA analysis to assess the reproducibility and the quality of the mRNA-seq dataset from Qin et al. (2022). Similarly, the biological repeats are not grouped according to the general biological expectation (Figure 3a). Mock libraries and Botrytis-infected libraries were grouped with each other (Figure 3a), indicating that the mRNA sequencing data have low repeatability and quality. To determine if the difference within each group is statistically significant, we performed differential expression analyses between infected and mock samples at each time point. We visualized results using volcano plots based on the FDR (false discovery rate) and foldchange value between different replicates and treatments. When using an FDR <0.05 and log2 fold change >1 as a cutoff, there are few (90 and 23) mRNAs with significant changes in "Infection 12 hours" and "Infection 24 hours," respectively ( Figure 3b-d). Thus, using these data to determine the expression levels of Bc-sRNAs targeted mRNAs is most likely inaccurate. The mRNA-seq datasets of ku70 and the ∆bcdcl1/∆bcdcl2/ku70 were still not available, therefore, we were not able to check whether the silencing effect of sRNA targets of tomato was eliminated in the dcl double mutant. First, estimating sRNA production per retrotransposon (RT) locus based on total sRNA read counts makes no sense because many sRNAs map multiple times and to different RT loci, leading to a doubling of sRNA read counts. Without knowing which element(s) is/are true sRNAproducing locus/i, the counts should be at least divided by the number of mapping events for estimation. Therefore, the 10% mentioned above of sRNA reads derived from the ms3003 locus is misleading.

sRNA production by the ms3003 element
Second, there are even more fragmented RT elements (loci) -we have previously identified additional 221 RT loci of > 400 bp fragments (Porquier et al., 2021) dispersed over the B. cinerea genome, that have not been considered as a source of sRNA production by Qin et al.,and which 13 are making a valid estimation of the relative sRNA production per individual RT element/locus extremely complicated (if not impossible).
Third, it is not indicated by Qin et al. (2022), which sRNA library was used for estimating the percentage of sRNAs produced per RT element (as given in table.S3). This is highly relevant because we recently found that sRNAs mapping to different RTs, in particular belonging to the BcGypsy3 family, are highly up-regulated during tomato infection (Porquier et al., 2021).
There are contradictory results and interpretations in Qin et al. (2022). Deleting the ms3003 element certainly did not lead to a drastic loss of sRNA production, as Qin et al. (2022) suggested.
In fact, their data support the continuous production of RT-derived sRNAs, as shown in figure 5c of Qin et al. In order to make any statement about altered sRNA production, the authors should provide sRNA sequencing data for the ms3003 mutants.

Reproducibility of infection assay with the ∆bcdcl1/∆bcdcl2/ku70 mutant strains
Qin et al. argue that the newly generated ∆bcdcl1/∆bcdcl2 mutant in the ku70 mutant background did not exhibit any reduction of virulence on tomato leaves and other tested plant species. For this experiment, they used a relatively high spore concentration of 1000 conidia/μl suspended in potato dextrose broth (PDB), a rich culture medium that allows rapid growth of B. cinerea in vitro, to inoculate the detached leaves. It is worth noting that tomato is a much more favored host for B. cinerea than Arabidopsis.
We have now repeated infection assays with detached tomato leaves using 10 μl of 10 conidia/μl spore concentrations, resuspended in 1% malt extract. We think this lower spore concentration better represents the natural infection pressure. Under this infection condition, we 14 observed significantly reduced lesion formation in two independent ∆bcdcl1/∆bcdcl2/ku70 mutant strains at 60 hours post-inoculation ( Figure 4A).
We would like to point out that the involvement of DCLs in fungal virulence was also observed in other fungal species and by other laboratories. Several studies reported that dcl deletion in fungal pathogens compromises infection, for example, in Fusarium graminearum, Penicillium italicum, Colletotrichum gloeosporioides and Valsa mali , Werner et al., 2021, Yin et al., 2020, Feng et al., 2017.

Qin et al. utilized an inappropriate NHEJ DNA repair mutant ku70 as a background to generate ∆bcdcl1/∆bcdcl2 mutants
Qin et al. claimed that they generated ∆bcdcl1/∆bcdcl2 mutant strains using the Cas9 RNP strategy (Leisen et al., 2022) in a B. cinerea B05.10 ku70 mutant background. This B. cinerea ku70 mutant strain is compromised in a non-homologous end-joining DNA repair (NHEJ) system and favors homologous recombination to gain targeted gene knock-out , Pinedo et al., 2008. However, double deletion of NHEJ and DCLs genes should be considered with caution for two reasons: First, a ku70 mutant may alter fungal pathogenicity. For example, the ku70 mutant of the fungus Penicillium digitatum exhibited altered growth and conidia production (Gandía et al., 2016), which means that pleiotropic effects may occur due to increased sensitivity to DNA damage.
In addition, RNAi mutants including DCLs showed increased sensitivity to DNA damage in Neurospora crassa (Lee et al., 2009), a phenomenon well-studied in the mammalian field as well (Tang & Ren, 2012, Bonath et al., 2018, Lu et al., 2018. Therefore, NHEJ plus DCLs deletion might likely cause unpredictable complications regarding DNA repair and genome integrity.

15
Indeed, PCA analysis showed that the sRNA profiles of the ku70 mutant had huge differences from that in wild type B. cinerea B05.10 strain ( Figure 4B).
Second, it might also be a problem when combining ku70 mutant with CRISPR CAS9 system. CAS9 is known to generate unpredicted off-target cleavages, which largely rely on NHEJ DNA repair system. The combination of CAS9 and ku70 mutation can generate more unpredicted mutation sites in the background. Many people in the field of mammalian and plant fungal pathogens are always hesitant in using ku70/ku80 mutants because they are fully aware that knockout NHEJ DNA repair function may cause many unpredicted mutations in the genome and also make the organisms more sensitive to DNA damage.  Table 1). Although accumulating in similar read counts, only Bc-siRNA3.1, BcsiR3.2 and Bc-siR5 bound to tomato AGO1, but not Bc-sRNA #6. Furthermore, the Bc-siRNA #3, which had relatively low read numbers, was still detected in tomato AGO1 (Figure 1b). These results highlight the importance of the tomato AGO1 co-IP sRNA analysis in order to identify potential cross-kingdom Bc-sRNAs.

DISCUSSION
Small RNA-mediated RNA interference (RNAi) is a conserved eukaryotic gene-silencing mechanism. Accumulating evidence has shown that sRNAs can not only regulate gene expression in the same cells, but they can also travel between interacting organisms to induce gene silencing in trans, which is called cross-kingdom or cross-species RNAi (Halder et al., 2022, Wong-Bajracharya et al., 2022, Cui et al., 2019, Ren et al., 2019, Dunker et al., 2020, Shahid et al., 2018, Ji et al., 2021, Zhang et al., 2022, Buck et al., 2014, Chow et al., 2019, Sahr et al., 2022. Recently, Qin et al. (2022)  cinerea samples. The distribution of sRNA lengths in the ku70 sample also differed from that of the cultured B05.10 B. cinerea, indicating a huge difference between these two strains.
Additionally, our experimental validation has confirmed that not all the transposon-derived sRNAs have been eliminated from the ∆bcdcl1/∆bcdcl2/ku70 triple mutant. Regarding the infection assay of ∆bcdcl1/∆bcdcl2/ku70, more conditions should be considered, for example, different spore concentration and different host species for infection need to be taken into consideration. We observed significantly reduced lesion formation in two independent ∆bcdcl1/∆bcdcl2/ku70 mutant strains using lower spore concentrations for infection assay at 60 hours post-inoculation.
In summary, cross-kingdom RNAi is present between B. cinerea and tomato. Some B. cinerea sRNAs can travel into tomato cells and bind to tomato AGO1 to silence target genes in the host.      Production of 9 BcsRNAs was tested in different B. cinerea genotypes by stem-loop RT-PCR.

Tomato AGO1 immunoprecipitation (IP)
Tomato AGO1 antibody from Dr. Borkovich's lab was used to pull down the sRNAs associated with the AGO1 protein (Ji et al., 2021). Tomato AGO1 IP was performed with 5 g fresh leaves collected at 0, 12, 16 and 24 h after spray inoculation with B. cinerea. Tissues were collected and ground in liquid nitrogen before resuspension in 5 ml of extraction buffer (20 mM Tris-HCl at pH 7.5, 150 mM NaCl, 5mM MgCl2, 1 mM DTT, 0.5% NP40, proteinase inhibitor cocktail; Sigma).
The extract was then precleared by incubation with 50 ul of Protein A-agarose beads (Roche) at 4°C for 30 min. The precleared extract was incubated with 10 ul of AGO1-specific antibody for 6 hours. A volume containing 50 ul of Protein A -agarose beads (Roche) was then added to the extract before overnight incubation. Immunoprecipitates were washed three times (5 min each) in wash buffer (20 mM Tris-HCl at pH 7.5, 150 mM NaCl, 5mM MgCl2, 1 mM DTT, 0.5% Triton x-100, proteinase inhibitor cocktail; Sigma). RNA was extracted from the bead-bound EVs with TRIzol reagent (Invitrogen). The sRNA RT-PCR was performed as previously described (Weiberg et al., 2013). Primer sequences are provided in Supplementary Table 2.

Bc-sRNA target gene expression analysis
The B. cinerea spores were diluted in 1% sabouraud maltose broth buffer to a final concentration of 1 × 10 5 spores/ml. The diluted spores were evenly sprayed on 5-to 6-week-old tomato plants.
The tomato leaves were collected after 0, 12, 16, and 24 hours of infection. Then total RNA was extracted from tomato leaf samples using TRIzol reagent (Invitrogen). The Bc-sRNA target gene expression was detected using quantitative real-time PCR. Primer sequences are provided in Supplementary Table 2.

Leaf inoculation
Tomato detached leaves were collected from 5-6 weeks old tomato plants. B. cinerea conidia were harvested from HA plate, where Botrytis had grown for 2 weeks. The conidia were adjusted the final concentration to 10 conidia/μl and resuspended in 1% malt extract for 1 hour. Tomato leaves were inoculated with 10 μl conidia suspension. Infected leaves were incubated in a humidity box.
Infected phenotypes were photographed and quantified the lesion size using the Fiji software, ImageJ version 2.1.0/1.53c (Schindelin et al., 2012).

24
Mapping and assigning to genome of origin Genome sequences for the Botrytis cinerea B05.10 (ASM83294v1) and Solanum lycopersicum SL3.0 where downloaded from Ensembl v55 (Yates et al., 2022), and unique sequences were mapped using bowtie1 (Langmead et al., 2009), allowing 0 or 1 mismatch, saving those that mapped or not to each genome but without further mapping information. Sequences were assigned to the genome with fewer mismatches, to "both" if the number of mismatches was the same and "none" if no hit was found. All sequences assigned only to B. cinerea were mapped again with bowtie1 using the parameters -v 1 -k 10000 --best --strata. The read counts assigned to each sequence in each sample were divided equally among all mapping locations.

Botrytis genome annotation
Genomic coordinates for protein-coding genes, tRNA, rRNA, other non-coding RNA classes and pseudogenes, were taken from the Ensembl GFF file. Transposable Element annotation was obtained by running RepeatModeler2 (Flynn et al., 2020) and RepeatMasker (Smit et al., 2013(Smit et al., -2015. To help assigning mapped sequences to unique annotations (and avoid overcounting), a non-redundant genomic annotation was built where all regions with overlapping annotations were reduced to a single annotation. The following preference was used: rRNA > tRNA > Transposable Elements > protein-coding exons > other ncRNA > introns > pseudogenes. The counts for all sequences in a sample (Figure 2c), or divided by length (Figure 2b), were assigned to the annotation category at the central base of each mapping location. Data processing and plotting was done in R (Team., 2022), making use of Biostrings (H. Pagès, 2021), GenomicRanges (Lawrence et al., 2013), Rsamtools (Morgan M, 2021) and rtracklayer (Lawrence et al., 2009).

RNA seq data analysis method 2-Jason Stajich
Bioinformatic analyses to examine read size distribution and genomic context were performed with scripts developed and available in this github repository https://github.com/stajichlab/Botrytis_sRNA_profile. Sequence reads were processed with fastp (Chen et al., 2018) to detect and trim adaptors and remove low quality reads. Trimmed sequences were aligned to Solanum lycopersicum ITAG4 (doi: 10.1101/767764v1) downloaded from SolNet solgenomics.net and Botrytis cinerea B05-10 genome downloaded from FungiDB (Basenko et al., 2018) using STAR aligner (Dobin et al., 2013). The resulting SAM alignments were converted to BAM files with samtools and processed with a custom python scripts 'process_BAM_sRNA.py' and 'get_BAM_sRNA_unique.py' which relied on Biopython (Cock et al., 2009) and pybedtools (Dale et al., 2011) generate size distribution and classifications of genomic features overlapped by the read alignments and visualized as plots with custom R scripts `plot_sRNA_sizes.R` and `plot_sRNA_sizes_uniq.R`.

RNA seq data analysis method 3-Qiang Cai
Data was acquired from the NCBI database (https://www.ncbi.nlm.nih.gov/) with accession number PRJNA496584. The data was detected quality by using FastQC v0.11.9. Adapter clipping, low-quality bases removal and reads size selection were applied to the data with Trim Galore v0.6.10. The clean data was used to align to the B. cinerea B05.10 genome ASM83294v1 assembly and S. lycopersicum genome SL3.0 assembly by using bowie v1.3.1. The reads unmapped to S. lycopersicum genome and exactly mapped to the B. cinerea B05.10 genome were used for subsequent analysis. After counting sRNA reads, the CPM (Counts Per Million mapped reads) was calculated by using R package edgeR v3.36.0. Principal component analysis (PCA) was applied for the above samples with the R package PCAtools v2.6.0 to obtain more insight into the 26 separation between samples of each treatment separately. The ggplot2 v 3.4.0 package was used to draw a relative image in R.

S. lycopersicum transcriptome analysis
Data was acquired from the NCBI database (https://www.ncbi.nlm.nih.gov/) with accession number PRJNA496584. The above sequencing reads were filtered using fastp v0.23.2 with adapter trimming and quality filtering function. The filtered reads were mapped to B. cinerea genome Principal component analysis (PCA) analysis was applied with the R package PCAtools v2.6.0 to reflect the repeatability between all samples. DESeq2 v1.34.0 was used to differentially express gene analysis and ggplot2 v3.4.0 was applied to draw the volcano plot. Through the above method, the differentially expressed genes which met the condition that the absolute value of log2fold change was more than 1 and the FDR value was less than 0.05 were identified.