Retrotransposons as pathogenicity factors of the plant pathogenic fungus Botrytis cinerea

Abstract

Background

Retrotransposons are genetic elements inducing mutations in all domains of life. Despite their detrimental effect, retrotransposons can become temporarily active during epigenetic reprogramming and cellular stress response, which may accelerate host genome evolution. In fungal pathogens, a positive role has been attributed to retrotransposons when shaping genome architecture and expression of genes encoding pathogenicity factors; thus, retrotransposons are known to influence pathogenicity.

Results

We uncover a hitherto unknown role of fungal retrotransposons as being pathogenicity factors, themselves. The aggressive fungal plant pathogen, Botrytis cinerea, is known to deliver some long-terminal repeat (LTR) deriving regulatory trans-species small RNAs (BcsRNAs) into plant cells to suppress host gene expression for infection. We find that naturally occurring, less aggressive B. cinerea strains possess considerably lower copy numbers of LTR retrotransposons and had lost retrotransposon BcsRNA production. Using a transgenic proof-of-concept approach, we reconstitute retrotransposon expression in a BcsRNA-lacking B. cinerea strain, which results in enhanced aggressiveness in a retrotransposon and BcsRNA expression-dependent manner. Moreover, retrotransposon expression in B. cinerea leads to suppression of plant defence-related genes during infection.

Conclusions

We propose that retrotransposons are pathogenicity factors that manipulate host plant gene expression by encoding trans-species BcsRNAs. Taken together, the novelty that retrotransposons are pathogenicity factors will have a broad impact on studies of host-microbe interactions and pathology.

Background

Retrotransposons are genetic elements inducing mutations in the bacteria, archaea and eukaryote domains of life [1, 2]. In eukaryotes, distinct classes of suppressive, cis-regulatory sRNAs, such as PIWI-associated piRNAs in Drosophila and nematodes [3] and heterochromatic small interfering RNAs in plants, are produced for retrotransposon control [4, 5]. Despite their detrimental effect, retrotransposons become temporarily active during epigenetic reprogramming [6] and cellular stress response [7, 8], which may accelerate host genome evolution [9,10,11,12]. In fungal pathogens, a positive role has been attributed to retrotransposons to shaping genome architecture and expression of genes encoding pathogenicity factors [13, 14]; thus, retrotransposons can influence pathogenicity.

B. cinerea can infect > 1400 plant species, including important crops such as tomato [15, 16]. For infection, B. cinerea small RNAs (BcsRNAs) translocate into plants and hijack the plant Argonaute (AGO)1/RNA-induced silencing complex (RISC) to suppress host immunity genes. Cross-kingdom sRNA effectors have also been reported in other fungal, oomycete, bacterial, and parasitic plant species [17,18,19,20], making cross-kingdom RNA interference (ckRNAi) a remarkably common phenomenon in plant-pathogen interactions [21]. The vast majority of BcsRNA effectors derive from long-terminal repeat (LTR) retrotransposons [22, 23], raising the hypothesis that retrotransposons could play a role as fungal pathogenicity factors.

Results

Botrytis cinerea strains carrying LTR retrotransposons are more aggressive

Previous transposon annotation in a published B. cinerea genome release (strain B05.10) revealed 83 full-length LTR retrotransposon copies separated into nine different consensus classes (Additional file 1: Figure S1) either belonging to the Gypsy or the Copia superfamily [24] (Additional file 2: Table S1). Using this B. cinerea LTR retrotransposon collection, we defined six subfamilies according to phylogenetic analysis, which we named BcGypsy1-BcGypsy4 and BcCopia1-BcCopia2 (Fig. 1a, Additional file 2: Tables S1, S2). We analysed previously published sRNA-seq data (raw data are available at NCBI GEO: GSE45323, GSE45321) regarding LTR retrotransposon BcsRNA production in B. cinerea axenic culture grown on agar plates or during infection of the host plants Solanum lycopersicum (tomato) and Arabidopsis thaliana [23]. LTR retrotransposon BcsRNA production in axenic culture indicated that only BcGypsy1, BcGypsy3 and BcGypsy4 produced significant amounts of BcsRNAs (Additional file 1: Figure S2a). During plant infection, more than 80% of LTR retrotransposon BcsRNAs mapped to BcGypsy3. These BcsRNAs displayed induced accumulation at early infection time point (Fig. 1b, Additional file 1: Figure S2b) in contrast to BcsRNAs derived from other genomic loci (Additional file 1: Figure S2c) suggesting a role in pathogenicity.

Fig. 1

Botrytis cinerea strains carrying LTR retrotransposons are more aggressive. a LTR retrotransposon phylogenetic relationship of six subfamilies, BcGypsy1-4, BcCopia1-2, identified in the strain B05.10. b LTR retrotransposon BcsRNA abundance of the six subfamilies (BcGypsy1-4, BcCopia1-2) at different time points of S. lycopersicum or A. thaliana host infection (raw data are available at NCBI GEO: GSE45323, GSE45321). c Growth phenotype of the six B. cinerea strains. Scale bars indicate 20 mm. d Genotyping PCR of the six LTR retrotransposon subfamilies in the six B. cinerea strains. e Pathogenicity assay of the six B. cinerea strains on tomato leaves quantifying microlesion area by Trypan Blue staining at 24 hpi. Microlesions in Trypan Blue images were quantified in eight leaf discs per strain as relative grey scale (rel. counts) in relation to the total leaf disc. Scale bars in leaf disc image indicate 1 mm and in higher magnification images 100 μm. f Pathogenicity assay of the six B. cinerea strains on tomato leaves indicating lesion area of > 20 infection sites at 48 hpi. Scale bar indicates 1 cm. For f, infection experiments were repeated three times with similar results. In e and f, a significant difference is indicated by letters and was tested by one-way ANOVA using Tukey HSD test with p < 0.05

To investigate the contribution of LTR retrotransposons to the pathogenicity of B. cinerea, we analysed six B. cinerea strains collected from different host plants and geographical origins (Additional file 2: Table S3) at genomic, sRNA transcriptomic, and disease phenotypic levels. Strains grown on nutrient-rich agar showed no obvious phenotypic differences; however, the strain D13_TF grew faster (Fig. 1c, Additional file 1: Figure S3a) and did not form sclerotia (Additional file 1: Figure S3b), a dormant tissue that enables B. cinerea to survive non-favourable conditions. We designed primers to genotype the six B. cinerea strains for the genomic presence of the six LTR retrotransposon subfamilies by PCR. The strains D08_H24, D14_KF and D13_TS were negatively tested for all LTR retrotransposon subfamilies. In contrast, B05.10, N11_KW and D13_TF were tested positive for all subfamilies (Fig. 1d). We performed pathogen assays with the six B. cinerea strains on tomato leaves to assess aggressiveness. We defined strain aggressiveness as a quantitative level of its ability to induce disease symptoms in the infected host. Tomato was chosen as a suitable host plant, as it was found to be susceptible to LTR retrotransposon BcsRNAs triggering ckRNAi [23]. The three strains negatively tested for LTR retrotransposons were less aggressive compared to the three LTR retrotransposon positively tested strains, considering induction of primary lesion formation at 24 h post inoculation (hpi) (Fig. 1e) and extended lesion area on infected leaves at 48 hpi (Fig. 1f). Based on these results, we found a positive relationship between LTR retrotransposons and B. cinerea aggressiveness.

Aggressive B. cinerea strains produce massive retrotransposon small RNAs

To gain further insights into the relationship between B. cinerea aggressiveness and LTR retrotransposon BcsRNAs, we performed a comparative sRNA-seq analysis using B. cinerea cultures. Remarkably, genome-wide BcsRNA mapping revealed that only the three most aggressive strains produced massive amounts of transposon-derived BcsRNAs (Fig. 2a). After filtering out rRNA reads, we further quantified relative read numbers of BcsRNAs aligning to annotated transfer (t)RNAs, small nuclear/nucleolar (sn/sno)RNAs, protein coding genes (mRNAs) or transposons in all six B. cinerea strains. Consistently, the three most aggressive strains produced high amounts of transposon-derived BcsRNAs between 2 and 16%, compared to 0.03–0.12% by the three less aggressive strains (Fig. 2b, Additional file 2: Table S4). In this analysis, transposon-derived BcsRNAs mostly mapped to LTR retrotransposons (> 99.4% in the most aggressive strains), leaving BcsRNAs that derived from DNA transposons (Additional file 2: Table S2) of minor proportion (Additional file 2: Table S4). Size profiles of total BcsRNAs from the six strains were rather diverse (Fig. 2c, Additional file 1: Figure S4a), but LTR retrotransposon BcsRNAs showed a clear size preference of 21-22 nucleotides (nt) and preferentially 5′ terminal Uracil (U) (Fig. 2d, Additional file 1: Figure S4a). In plants, AGO1 commonly associates with endogenous 21-22 nt and 5′U sRNAs to form an AGO1/RISC [25], which likely explains why B. cinerea LTR retrotransposon BcsRNAs bind to the plant AGO1 and induce ckRNAi [22, 23]. Most LTR retrotransposon BcsRNAs were found to be unique or were share by two of the more aggressive B. cinerea strains (Additional file 1: Figure S4b). In this regard, the majority of LTR retrotransposon BcsRNA effectors (29 out of 48) mapped to BcGypsy1 and BcGypsy3 subfamilies. In accordance, the less aggressive strains D08_H24, D14_KF and D13_TS showed nearly no expressed BcsRNAs aligning to BcGypsy1 and BcGypsy3, whereas many were found in the aggressive B05.10, N11_KW and D13_TF strains (Additional file 1: Figure S5). The strain D13_TF produced less LTR retrotransposon BcsRNAs compared to B05.10 and N11_KW (Fig. 2a, b, Additional file 1: Figure S5); however, it induced in average the biggest lesions (Fig. 1f). The relative large lesion formation could be explained with its fast growth (Fig. 1c, Additional file 1: Figure S3a). With these results, we found evidence that the less aggressive B. cinerea strains lacked production of LTR retrotransposon BcsRNAs, which further supported that LTR retrotransposon BcsRNAs are important for pathogenicity.

Fig. 2

Comparative sRNA-seq analysis of the six B. cinerea strains. a Genome-wide BcsRNA maps of the six B. cinerea strains, coverage represented as log(RPM). b Relative composition of BcsRNAs mapping to distinct genomic loci with numbers giving the relative percentage. c Size profiles of total BcsRNAs in the six B. cinerea strains. d Sizes profiles and 5′ first nucleotide distribution of LTR retrotransposon BcsRNAs in the six B. cinerea strains; in c and d reads per million (RPM)

Less aggressive B. cinerea strains carry only silenced LTR retrotransposon relics

To clarify if the reason for the absence of LTR retrotransposon BcsRNAs in less aggressive B. cinerea strains was due to the lack of LTR retrotransposons, as indicated by genotyping PCR (Fig. 1d), we re-sequenced the genomes of the two less aggressive strains D08_H24 and D14_FK using hybrid Nanopore long-read in combination with Illumina short-read sequencing. Scaffold assembly resulted in 49 contigs for D08_H24 and 28 contigs for D14_KF. Pairwise genome alignment and synteny analysis of D08_H24 and D14_KF with the published B05.10 genome [26] revealed coverage of all 18 chromosomes (Fig. 3a). The nearly gapless full genome assemblies of D08_H24 and D14_KF allowed us to perform transposon annotation and analysis. Using the REPET pipeline [27] on the D08_H24 and D14_KF genomes, we did find no Gypsy LTR retrotransposon and only one Copia LTR retrotransposon with 14 and 10 full-length copies, respectively (Fig. 3b, Additional file 2: Table S5). To complement our REPET analysis, we performed Blastn search in the published B05.10, and the D08_H24 and D14_KF re-sequenced genomes, using REPET-annotated full-length LTR retrotransposons of B05.10 as queries and allowing a minimum alignment length of 400 nucleotides (nt). By this step, we identified further Gypsy and Copia elements in B05.10, D08_H24 and D14_KF (Fig. 3b), with fewer and more truncated ones in D08_H24 and D14_KF compared to in B05.10 (Fig. 3c, Additional file 1: Figure S6a-b). Using all retrotransposon copies annotated in B05.10 and the re-sequenced D08_H24 and D14_KF genomes either identified by REPET or Blastn for sRNA-seq read mapping, we confirmed that D08_H24 and D14_KF produced marginal amounts of LTR retrotransposon BcsRNAs compared to B05.10 (Fig. 3b, Additional file 2: Table S6). Since the reduced accumulation of retrotransposons in the strains D08_H24 and D14_KF was not reflected in their genome sizes (Fig. 3b), we would presume that other types of DNA have compensated for the reduction of transposons. A striking difference regarding LTR retrotransposons between D08_H24, D14_KF, and B05.10 was the GC content. While D08_H24 and D14_KF carried exclusively full-length retrotransposons with low GC content (< 30%), B05.10 displayed two distinct fractions of low GC (< 30%) or high GC (> 40%) contents (Fig. 3d, Additional file 2: Tables S2, S5). We obtained similar results for truncated retrotransposon copies identified by Blastn, although few truncated copies in D08_H24 and D14_KF showed GC content > 40% (Additional file 1: Figure S7). Low GC content let us presume that LTR retrotransposons were mutated by the Repeat-induced point mutation (RIP), which is a fungal genome surveillance mechanism mutating C to T bases within duplicated sequences [28]. Previous RIP analysis in B. cinerea had indicated bias for CA to TA and CT to TT dinucleotides [29]. In order to estimate the occurrence of RIP on REPET-detected full-length retrotransposons, we calculated the TA/AT index. We found diverse TA/AT ratios in B05.10, but only high (> 0.89) TA/AT index in D08_H24 and D14_KF (Fig. 3e, Additional file 2: Table S7), a threshold considered to indicate RIP [30]. Similarly, comparative analysis of both full-length and truncated retrotransposon copies in B05.10, D08_H24 and D14_KF using the RIPCAL program [30] revealed RIP in all copies of D08_H24 and D14_KF, but some copies in B05.10 with low RIP, as exemplified for the BcGypsy3 subfamily (Additional file 1: Figure S8). All full-length LTR retrotransposons found in D14_KF or D08_H24 displayed index and content (< 30%). This finding possibly explains why LTR retrotransposon genotyping PCR was negative for D08_H24 and D14_KF (Fig. 1d), as primer sequences did not match the mutated LTR retrotransposon sequences. In fungi, RIP often leads to transcriptional silencing [31], and we anticipated that loss of LTR retrotransposon BcsRNAs could be a result of low LTR retrotransposon expression. Indeed, expression values of BcsRNAs mapping to BcGypsy1, BcGypsy3 and BcGypsy4 correlated positively with the LTR retrotransposon GC content and negatively with their TA/AT index (Fig. 3f, Additional file 2: Table S8) in B05.10. GC-rich copies of BcGypsy2, BcCopia1 and BcCopia2 showed low BcsRNA mapping accumulation in general; still, more reads mapped to GC-rich copies within the corresponding LTR retrotransposon subfamily (Additional file 2: Table S9). Our genome-wide comparative analysis of LTR retrotransposons revealed that less aggressive B. cinerea strains possessed fewer copy numbers with low GC content, which did not produce LTR retrotransposon BcsRNAs.

Fig. 3

Comparative genome analysis of LTR retrotransposons. a Whole genome chromosomal alignment analysis. b–e Comparative analysis of LTR retrotransposons identified in the strains B05.10, D08_H24 and D14_KF, and differences in the level of copy numbers, truncation, GC content (%), and the TA/AT dinucleotide ratio with threshold line shown at 0.89. f Correlation analysis between LTR retrotransposon BcsRNA read abundance and GC content (%) or TA/AT dinucleotide ratio of the strain B05.10. r gives the Pearson correlation coefficient

The B. cinerea retrotransposon BcGypsy3 promotes host plant infection

Having observed a positive relationship between full-length, GC-rich LTR retrotransposons, LTR retrotransposon BcsRNAs, and pathogenicity, we sought to validate the concept of retrotransposon as pathogenicity factor by a transgenic approach. We cloned a GC-rich BcGypsy3 from the aggressive B05.10 strain to transform the less aggressive D08_H24, anticipating that transformed D08_H24 would produce BcGypsy3 BcsRNAs and gain pathogenicity. We chose a GC-rich BcGypsy3 as a major source of BcsRNAs, and chose D08_H24, because this strain did not possess any full-length BcGypsy3 (Additional file 1: Figures S6b, S9) and did not produce significant amounts of LTR retrotransposon BcsRNAs (Fig. 2, Additional file 2: Table S6). A 4.98 kilobase BcGypsy3 was cloned into a fungal expression vector. We excluded flanking LTR DNA from cloning, because these typically contain binding sites of tRNA primers that initiate reverse transcription and transposon transposition [32], and because BcsRNAs did not derive from the LTR regions (Additional file 1: Figure S2b). We generated two BcGypsy3 cassettes driven by constitutive promoters to either produce single strand (ss)BcGypsy3 sense RNAs or sense and antisense RNAs to form double strand (ds)BcGypsy3 RNAs to enhance BcsRNA production. Empty vector (EV) was used as transformation control (Additional file 1: Figure S10a). Transgenes were inserted into the Nitrate reductase D (BcniaD) locus by homologous recombination, which was previously established for targeted transgene insertion in B. cinerea without affecting pathogenicity [33]. Genomic insertion and BcGypsy3 mRNA expression were validated by PCR (Additional file 1: Figure S10b) and quantitative reverse transcription (qRT)-PCR (Fig. 4a), respectively. We performed pathogen assays with BcGypsy3 transformants that grew comparably on nutrient-rich agar (Additional file 1: Figure S10c) to test for increased aggressiveness. Four independent dsBcGypsy3 transformants induced significantly larger lesions than an EV transformant (Fig. 4b, Additional file 1: Figure S11). Replicating this experiment with four ssBcGypsy3 transformants gave similar results, albeit gained pathogenicity was not as strong as in the dsBcGypsy3 transformants (Fig. 4b, Additional file 1: Figure S11). These results were in line with a higher BcGypsy3 mRNA expression found in the dsBcGypsy3 transformants compared to the ssBcGypsy3 transformants (Fig. 4a). We repeated this experiment with two isolated BcGypsy3 transformants carrying a 258 bp 5′ LTR element instead of a promoter to simulate a more natural situation and obtained again larger lesion formation compared to the EV control strain (Additional file 1: Figure S12). We selected one ssBcGypsy3 (transformant #51), one dsBcGypsy3 (transformant #56) and the EV strain (transformant #18) for sRNA-seq analysis to relate BcGypsy3 BcsRNA production to pathogenicity. BcGypsy3 BcsRNA accumulation was moderately higher in the dsBcGypsy3 compared to the ssBcGypsy3 (Fig. 4d), which was in line with the higher BcGypsy3 mRNA expression and aggressiveness of dsBcGypsy3 transformants. However, overall BcGypsy3 BcsRNA accumulation in the transgenic D08_H24 was not as high as in B05.10, although the dsBcGypsy3 transformants displayed similar BcGypsy3 mRNA expression (Fig. 4a). We therefore suspected that insertion of the BcGypsy3 transgene into the BcniaD locus might have limited BcsRNA production. To test this possibility, we repeated transformation without BcniaD homologous recombination flanking DNA for BcGypsy3 random genome integration. We isolated two independent, random-inserted BcGypsy3 (randBcGypsy3) transformants #1, #3 with confirmed BcGypsy3 mRNA expression (Fig. 4a) and exhibiting comparable growth to the other isolated transformants (Additional file 1: Figure S10c). The randBcGypsy3 #1 was significantly more aggressive compared to D08_H24 wild type (wt) and displayed the highest BcGypsy3 mRNA expression (Fig. 4a, c). Sequencing sRNAs of the randBcGypsy3#1 revealed that production of BcGypsy3 BcsRNAs was at approximately 10 times higher and more diverse than in the transformants with BcGypsy3 inserted into the BcniaD locus (Fig. 4d, Additional file 1: Figure S13), which correlated with the highest BcGypsy3 mRNA expression in the randBcGypsy3 #1 (Fig. 4a). However, BcGypsy3 BcsRNA production in the randBcGypsy3 #1 was around 35 times lower than in B05.10 (Fig. 4d). Taken together, transgene-induced BcGypsy3 BcsRNA expression led to enhanced aggressiveness, which strongly supported that this retrotransposon is a pathogenicity factor in B. cinerea.

Fig. 4

BcGypsy3 promotes host plant infection. a BcGypsy3 mRNA expression in single strand (ss)BcGypsy3 (#4, #49, #50, #51), double strand (ds)BcGypsy3 (#5, #6, #56, #57) or empty vector (EV) (#18) transformants with the BcGypsy3 transgene inserted in the BcniaD locus of D08_H24. Further, randomly inserted ssBcGypsy3 (randBcGypsy3) tranformants #1, #3 into the D08_H24 genome, as well as the wild type strains D08_H24 and B05.10 are shown. Data points represent three biological replicates. A significant difference as indicated by letters was tested by one-way ANOVA using Tukey HSD test with p < 0.05. b Pathogenicity assay with dropped spore suspension of BcGypsy3 transformants on tomato leaves quantifying lesion area of > 20 infection sites at 48 hpi. A significant difference as indicated by letters was tested by one-way ANOVA using Tukey HSD test with p < 0.05. c Pathogenicity assay with agar plugs of randBcGypsy3 transformants on tomato leaves quantifying lesion area of minimum 8 infection sites at 48 hpi. Significant difference was tested by two-sided Student’s t test. For b and c, similar results were obtained in two independent infection experiments. d BcsRNA-seq analysis at the BcGypsy3 transgene locus representing ssBcGypsy3 #51, dsBcGypsy3 #56, randBcGypsy3 #1 or EV #18 transformants with blue bars indicate sense and red bars antisense reads. Table shows normalized read counts (RPM) of BcGypsy3 BcsRNAs. e On top, stem-loop RT-PCR showed BcsRNA3.1, BcsRNA3.2 and BcsRNA20 expression in the randBcGypsy3 #1. BcGypsy3 mRNA expression and BcTub were used as controls. On the bottom, qRT-PCR of S. lycopersicum target mRNA expression after water treatment (mock), infection with D08_H24 wt or randBcGypsy3 #1. Each data point represents a biological replicate. f On the top, scheme of RNA-seq experiment infecting A. thaliana either with D08_H24 wt or randBcGypsy3 #1 comprising 12 biological replicates for each treatment. On the bottom, heat map showing differentially expressed A. thaliana genes comparing infection with D08_H24 wt or randBcGypsy3 #1

We hypothesized that the BcGypsy3 transformants exhibited enhanced aggressiveness due to the production of BcGypsy3 BcsRNAs suppressing plant target genes. To test this hypothesis, we chose three known BcGypsy3 BcsRNA effectors, BcsRNA3.1, BcsRNA3.2 and BcsRNA20, and their target genes in tomato (Additional file 1: Figure S14a) [23]. We confirmed expression of BcsRNA3.1, BcsRNA3.2 and BcsRNA20 in the randBcGypsy3 #1 by stem-loop RT-PCR (Fig. 4e). We then infected plants either with randBcGypsy3 #1 or D08_H24 wt in comparison to mock-treated plants and measured mRNA levels of the tomato Mitogen-activated protein kinase kinase kinase 4 (SlMPKKK4), Class E vacuolar protein-sorting machinery protein hse1 (Slhse1, Solyc09g014790) and the Basic helix-loop-helix (Bhlh)63 transcription factor (Solyc03g120530). The tomato target genes showed downregulation upon infection with randBcGypsy3 #1 compared to D08_H24 wt (Fig. 4e). To repeat this experiment with a different plant species, we infected A. thaliana with D08_H24 wt or randBcGypsy3 #1 to quantify gene expression of the known BcsRNA3.1 and BcsRNA3.2 target genes AtMPK1, AtMPK2 and a Cell wall-associated kinase (AtWAK) (Additional file 1: Figure S14b) [23]. AtMPK1 showed significant down-regulation, and AtMPK2 in tendency, upon infection with randBcGypsy3 #1 compared to D08_H24 wt (Additional file 1: Figure S14c). The A. thaliana infection-responsive genes AtPlant Defensin (PDF)1.2 and AtPathogensis-related protein (PR)1, which were all not predicted targets of BcsRNAs, did not show downregulation in plants infected with randBcGypsy3 #1 compared to D08_H24 wt (Additional file 1: Figure S14c) indicating that plant gene downregulation was not a general effect of randBcGypsy3 infection. The downregulation of target genes supported the role of BcGypsy3 BcsRNAs in plant gene manipulation.

To further explore the impact of the BcGypsy3 transgene on the plant mRNA transcriptome during infection, we conducted an RNA-seq experiment comparing A. thaliana plants infected with randBcGypsy3 #1 or D08_H24 wt (Fig. 4f). Differential gene expression analysis of A. thaliana mRNAs indicated only moderate effects on a small subset of A. thaliana genes. In total, we identified 7 up- and 8 downregulated candidate genes when plants were infected with the randBcGypsy3 #1 by differential expression analysis applying a false discovery rate (FDR) cut-off at 0.05 (Fig. 4f, Additional file 2: Table S10). In consistence, a similar pattern of up- or downregulation was evident for the candidate genes according to the Bio-Analytic Resource for Plant Biology (BAR) database [34] (Additional file 1: Figure S15). Down-regulated genes included the Small auxin upregulated RNA 78 (AT1G72430) and the Catalase (CAT)2 (AT4G35090), both related to auxin signalling. Repression of auxin response was previously shown to increase A. thaliana susceptibility to B. cinerea [35]. CAT2 as well as the downregulated cysteine protease Response to dehydration 21 (AT1G47128) are A. thaliana resistance factors against B. cinerea [36, 37]. Upregulated genes related to stress response included the transcriptional repressor NF-X-LIKE 1 (At1G10170) and the transcriptional repressor Jasmonate-Zim-Domain (JAZ)5 (At1G17380), both being negative regulators of A. thaliana defence response against fungal phytotoxins [38] and B. cinerea [39], respectively. Thus, BcGypsy3 contributed to manipulating plant gene expression, possibly via BcGypsy3 BcsRNA-induced ckRNAi. In this sense, we predicted 5 out of the 7 randBcGypsy3 #1 downregulated A. thaliana genes as targets of BcGypsy3 BcsRNAs by psRNATarget [40], but also found predicted BcsRNA targets among the up-regulated A. thaliana genes suggesting both direct and indirect effects on plant gene regulation (Additional file 2: Table S11).

Discussion

Our data provide lines of evidence that retrotransposons are pathogenicity factors in B. cinerea. Previous observation suggested that B. cinerea strains containing transposons colonize more frequently host plants and might induce stronger disease severity compared to transposon-free strains [41, 42]. In this regard, we previously discovered that B. cinerea delivers BcsRNAs into plants to suppress host immunity genes [23, 43] most of them derive from LTR retrotransposons [22]. We here showed that natural B. cinerea strains that have lost production of LTR retrotransposon BcsRNAs are less aggressive on tomato. We found that less aggressive B. cinerea strains failed to produce LTR retrotransposon BcsRNAs possibly due to RIP. RIP on LTR retrotransposons is common in ascomycete species [31]; however, it remains unclear why some LTR retrotransposons escaped RIP in the aggressive B. cinerea isolate B05.10. This could be due to positive selection inferred by host plant infection. Moreover, transgene-induced production of retrotransposon BcsRNAs in a non-retrotransposon BcsRNAs producing B. cinerea strain led to enhanced aggressiveness in a BcGypsy3 mRNA and BcGypsy3 BcsRNA expression-dependent manner. The level of transgene-induced BcGypsy3 BcsRNA production in the recipient strain D08_H24 was dependent on the genome insertion site, but did not fully restore BcsRNA production as found in the donor strain B05.10, although BcGypsy3 mRNA expression levels were similar. Unexpectedly, there was a positive correlation between BcGypsy3 BcsRNA expression and BcGypsy3 mRNA transcript levels. We speculate that post-transcriptional regulation might be the reason for moderate BcsRNA production, or that stronger BcsRNA production only occurs after several generations upon transformation [32, 44]. In addition, the retrotransposon transgene was transformed into D08_H24 without LTRs, which could be a hint that BcsRNA production also depends on LTRs. Performing RNA-seq with A. thaliana infected with D08_H24 wt or a BcGypsy3-transformant, we revealed first insights into the regulatory effects of B. cinerea retrotransposon BcsRNAs on the host plant transcriptome, comprising upregulation of immunity repressors as well as downregulation of genes involved in activating immunity against B. cinerea. Since this analysis gave only a snapshot on plant gene manipulation induced by a single retrotransposon, we predict that more and other plant genes might be regulated through diverse LTR retrotransposon BcsRNAs produced by the different B. cinerea isolates. Upregulation of LTR retrotransposon sRNAs was also evident in Magnaporthe oryzae, the fungal plant pathogen causing the rice blast disease, during infection of rice plants, as well as de-repression of LTR retrotransposon transcription was observed in the fungal wheat pathogen Zymoseptoria tritici [13] during wheat infection, suggesting a positive role of LTR retrotransposons in pathogenicity for other fungal pathogens. Moreover, animal-parasitic nematodes, such as Heligmosomoides polygyrus, secrete retrotransposon-derived sRNAs in extracellular vesicles [45], which are internalized by host cells and probably induce suppression of host mRNAs [46].

Conclusion

In summary, we found evidence for a correlation between pathogenicity and retrotransposons in a collection of natural fungal pathogen isolates. By a transgenic approach, we proved that retrotransposon expression in a weakly aggressive pathogen isolate significantly enhanced its aggressiveness level and that retrotransposons contribute to pathogenicity by producing trans-species sRNAs that manipulate host plant gene expression. It will be interesting to explore how common diverse pathogens and parasites utilize transposons as pathogenicity factors to colonize their hosts.

Methods

Strains materials, growth media and condition

The Botrytis cinerea (Pers.: Fr.) strain B05.10 as well as the five wild isolates D08_H24, D14_KF, D13_TF, D13_TS and N11_KW were used for this study. Routine cultivation was carried out on rich medium (RM; 10 g/L malt extract, 4 g/L yeast extract, 4 g/L glucose, 15 g/L agar) supplemented for mutant strains with hygromycin B (Sigma-Aldrich; 70 μg/mL) or nourseothricin (Werner Bioagents; 120 μg/mL). The plates were incubated on the bench in transparent plastic boxes under natural light conditions. Arabidopsis thaliana (L.) were grown on soil under short day condition (8 h light/ 16 h dark, 22 °C, 60% relative humidity). Solanum lycopersicum (L.) (tomato) cultivars Moneymaker or Heinz were grown under 16 h light/8 h dark, 24 °C, 60% humidity condition.

Pathogenicity assay

Pathogenicity assays were performed on detached tomato leaves from 4- to 5-week-old plants. B. cinerea conidia were resuspended in 1% malt extract at a final concentration of 5 × 10⁴ conidia/mL. Tomato leaves were inoculated with 20 μl conidia solution or agar plugs from 4-day-old Botrytis mycelia grown on RM agar, and inoculated leaves were kept in a humidity box. Infected leaves were photographed and lesion area was measured using the Fiji software (ImageJ version 2.1.0/1.53c).

Trypan Blue staining and microscopy

Infected tomato leaves were stained with Trypan Blue as described previously [47]. Microscopic images were taken with a DFC450 CCD-Camera (Leica) on a CTR 6000 microscope (Leica Microsystems). Microlesion quantification was done by the FIJI software by choosing the image-type 8-bit format (grey scale), scaling whole leaf diameter to 1, auto-setting of image adjust threshold and measure pixel counts in lesion area versus whole leaf area.

Genotyping PCR

Genomic DNA was isolated using the CTAB method followed by chloroform extraction and isopropanol precipitation [48]. PCR primers used are given in the Additional file 2: Table S12.

RNA extraction, Reverse transcription and quantitative PCR

Total RNA was isolated using a CTAB-based method [49]. Genomic DNA was removed using DNase I (Sigma-Aldrich) treatment and cDNA synthesis was performed with 1 μg total RNA using SuperScriptIII RT (ThermoFisher Scientific). Gene expression was measured by qPCR on a Quantstudio5 cycler (ThermoFischer Scientific) using the Primaquant low ROX qPCR master mix (Steinbrenner Laborsysteme). Primers are listed in the Additional file 2: Table S12. BcTubulin was used for normalization. Differential expression was calculated using the 2^-ΔΔCt method [50].

Stem-loop RT-PCR

Detection of BcsRNAs was carried out following the stem-loop RT-PCR protocol [51] using 1 μg of total RNA. Primers are listed in the Additional file 2: Table S12. The stem-loop RT-PCR products were visualized on 10% non-denaturing polyacrylamide gels.

Cloning and B. cinerea transformation

The different cloning constructs used for B. cinerea transformation were generated using the Golden Gate strategy [52] and are presented in the Additional file 1: Figure S9a, and used primers are listed in the Additional file 2: Table S12. Before fungal transformation, the plasmids were digested with the restriction endonuclease BsaI (NEB Biolabs) to isolate the transformation cassette of interest. The transformation of B. cinerea strain D08_H24 was performed as described before [53] with the following minor modifications. Transformed protoplasts were mixed into SH agar without antibiotics and incubated in dark for 24 hours. Upon pre-incubation, SH agar with protoplasts was covered with fresh SH agar containing 70 μg/mL of hygromycin B (Sigma-Aldrich) or 120 μg/ml of nourseothricin (Werner Bioagents) and the plates were further incubated in the dark until the isolation of putative mutants. Transformants with targeted insertion of transgenes into the BcniaD locus were further selected by plating on Czapek-Dox agar containing 0.4 M potassium chlorate as previously described [33]. Successful target insertion was confirmed by genomic PCR analysis.

Illumina sRNA-seq and data analysis

sRNAs were isolated from B. cinerea grown on RM for 4 days for high throughput sequencing, as previously described [23]. sRNAs were cloned for Illumina sequencing using the Next® Small RNA Prep kit (NEB) and sequenced on an Illumina HiSeq1500 platform. The Illumina sequencing data were analysed using the GALAXY Biostar server [54]. Raw data were de-multiplexed (Illumina Demultiplex, Galaxy Version 1.0.0) [40] and adapter sequences were removed (Clip adaptor sequence, Galaxy Version 1.0.0). Sequence raw data are deposited at the NCBI SRA server (BioProject PRJNA730711, SRA accessions SRR14576251- SRR14576256 and SRR14576386- SRR14576389). Reads were mapped to B. cinerea B05.10 reference genome assembly (BioProject: PRJNA15632) or de-novo genome assemblies of the strains D08_H24 and D14_KF using the BOWTIE algorithm (Galaxy Version 1.1.0) allowing zero mismatches (-v 0). Ribosomal RNAs (rRNA) were filtered out using the BOWTIE algorithm allowing three mismatches (-v 3). Reads were sorted to transfer RNA (tRNA), small nuclear/nucleolar RNA (snoRNA), messenger RNA, and transposon RNA reads using BOWTIE2 [55] with default settings. Upon classification, reads were count and normalized on total B. cinerea reads per million (RPM). Mapping files were visualized using a custom script based on ggplot2 (version 3.2.1) package in R (version 3.6.1). Target gene prediction of sRNAs was performed with the TAPIR program using a maximal score of 4.5 and a free energy ratio of 0.7 as thresholds [56] and psRNATarget with default settings [40].

RNA-seq

Conidia from 10-day-old mycelium cultures were harvested and diluted to 10⁶ conidia/mL in 1% malt extract (pH 5.7). Six-week-old Arabidopsis thaliana Col-0 plants were infected with those conidia (24 plants/strain dispatched into 2 humidity boxes). For each plant, 5 leaves were infected with 4 drops of 5 μl. After 24 h of incubation on the bench, 3 leaves from each plant were randomly harvested. Leaves containing the 4 spots of infection were cut and the infections from 2 leaves were pooled into a single Eppendorf tube and snap frozen in liquid nitrogen constituting one biological replicate. Total RNAs were extracted using a CTAB-based method [49]. Total RNA was subjected to mRNA sequencing using a version of the prime-seq method online available at https://www.protocols.io/view/prime-seq-s9veh66. This protocol is based on single-cell RNA-seq [57] and is a three-prime counting method that includes a sample-specific barcode sequence and unique molecular identifiers (UMI) for accurate quantification of gene expression. Illumina paired-end sequencing was performed on an HiSeq 1500 instrument. Raw data was demultiplexed using deML [58], adapters and poly-A tails were trimmed using cutadapt (version 2.3) and further preprocessed using the zUMIs pipeline [59] with STAR [60]. Processed raw data are deposited at the NCBI SRA server (BioProject PRJNA730711, SRA accessions SRR14577137-SRR14577147 and SRR14590656-SRR14590667). Reads were mapped to combined A. thaliana (TAIR10) or B. cinerea (ASM83294v1) genomes with Araport11 and ASM83294v1.41 gene annotations. Differential gene expression analysis was done in iDEP.91 [61] with FDR cut-off at 0.05 using limma-voom normalization.

Whole genome sequencing

WGS of B. cinerea D08_H24 and D14_KF strains was accomplished by hybrid Oxford Nanopore (Promethion) long-read sequencing and Illumina HiSeq1500 short-read sequencing. After base calling with guppy (version 2.3.7) and adapter removal with porechop (https://github.com/rrwick/Porechop), the long-reads where assembled using wtdbg2 [62] and polished with two rounds of racon [63]. The short reads were used to polish the assembly with pilon [64]. The Circos plotting library [65] was used via R circlize [66] to visualize whole genome comparisons. Assembled genome sequence contigs of D08_H24 and D14_KF strains are deposited at the NCBI SRA server (BioProject PRJNA730711).

Transposon annotation and RIP analysis

Transposable elements in the B. cinerea strain B05.10 were previously annotated [24] and used in this study. Annotation of transposable elements in D08_H24 and D14_KF genomes was performed using the REPET package as described previously [24]. Sequence alignment and phylogenetic tree generation was done using the CLC main workbench software (version 20.0.4, Qiagen) using the Neighbor Joining method and Jukes-Cantor nucleotide distance measure. Bootstrap analysis was run with 500 replicates. Blastn search was used to identify truncated retrotransposons (> 400 bp) with the following optional parameters -word_size 20 -max_target_seqs 50000 -gapopen 5 -gapextend 2 -reward 2 -penalty -3 -dust no -soft_masking false. The GC richest copy of each LTR retrotransposon subfamilies in B05.10 was chosen as a query. The sequences of the hits longer than 1000 bp were extracted and aligned with the refalign tool from the REPET package (https://urgi.versailles.inra.fr/Tools/REPET) to further run RIPCAL [30] analysis using the GC-richest query copy as a reference. RIPCAL was used to calculate all possible RIP mutations (CA ➔ TA, CT ➔ TT, CC ➔ TC, CG ➔ TG and their reverse complements). For the calculation of the RIP index TA/AT [67], only full-length copies from the REPET analysis were used. The different dinucleotide combinations were counted with the compseq tool (from the EMBOSS package) and used to calculate the RIP index.

Prediction of protein domains in LTR retrotransposons was done on a GC-rich copy sequence of each consensus class using the NCBI conserved domain database [68]. Long terminal repeats sequence and size were obtained from previous analysis [24]. For BcGypsy1, the additional ORF sequence (brtn) was obtained from Zhao et al. [69]. It is to note that this gene is not found in the last annotation of B05.10 genome. For the P26.1 consensus, the beginning of brtn is not present as in the other consensus grouped as the BcGypsy1 subfamily but a putative ORF is still present (dashed arrow in Fig.S1). The other represented additional ORFs were predicted in the last annotation of B05.10 genome (http://fungi.ensembl.org/Botrytis_cinerea). The predicted protein domain figure was generated using IBS (http://ibs.biocuckoo.org/).

Plots and statistical analysis

Plots and statistical tests were carried out using R studio (version 1.0.136, rstudio.com), ggplot2 (version 3.2.1) or Excel.