# Characterization of siRNAs clusters in Arabidopsis thaliana galls induced by the root-knot nematode Meloidogyne incognita

Clémence Medina, Martine da Rocha, Marc Magliano, Alizée Raptopoulo, Nathalie Marteu, Kevin Lebrigand, Pierre Abad, Bruno Favery, Stéphanie Jaubert-Possamai

► To cite this version:

Clémence Medina, Martine da Rocha, Marc Magliano, Alizée Raptopoulo, Nathalie Marteu, et al. Characterization of siRNAs clusters in Arabidopsis thaliana galls induced by the root-knot nematode Meloidogyne incognita. BMC Genomics, 2018, 19, 10.1186/s12864-018-5296-3. hal-02628557

HAL Id: hal-02628557

https://hal.inrae.fr/hal-02628557

Submitted on 26 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

Open Access

<figure><img src="image_2.png" /></figure>

CrossMark

# Characterization of siRNAs clusters in *Arabidopsis thaliana* galls induced by the root-knot nematode *Meloidogyne incognita*

Clémence Medina¹†, Martine da Rocha¹†, Marc Magliano¹, Alizée Raptopoulo¹, Nathalie Marteu¹, Kevin Lebrigand², Pierre Abad¹, Bruno Favery¹ and Stéphanie Jaubert-Possamai¹* ID

**Abstract**

**Background:** Root-knot nematodes (RKN), genus *Meloidogyne*, are plant parasitic worms that have the ability to transform root vascular cylinder cells into hypertrophied, multinucleate and metabolically over-active feeding cells. Redifferentiation into feeding cells is the result of a massive transcriptional reprogramming of root cells targeted by RKN. Since RKN are able to induce similar feeding cells in roots of thousands of plant species, these worms are thought to manipulate essential and conserved plant molecular pathways.

**Results:** Small non-coding RNAs of uninfected roots and infected root galls induced by *M. incognita* from *Arabidopsis thaliana* were sequenced by high throughput sequencing. siRNA populations were analysed by using the Shortstack algorithm. We identified siRNA clusters that are differentially expressed in infected roots and evidenced an over-representation of the 23–24 nt siRNAs in infected tissue. This size corresponds to heterochromatic siRNAs (hc-siRNAs) which are known to regulate expression of transposons and genes at the transcriptional level, mainly by inducing DNA methylation.

**Conclusions:** Correlation of siRNA clusters expression profile with transcriptomic data identified several protein coding genes that are candidates to be regulated by siRNAs at the transcriptional level by RNA directed DNA methylation (RdDM) pathway either directly or indirectly via silencing of neighbouring transposable elements.

**Keywords:** Gall, Giant cell, siRNA, Plant parasitic nematodes, Small RNA, Transcriptome regulation, Transposable element

**Background**

Small RNAs (sRNAs) are 20- to 24-nucleotide (nt) non-coding RNAs that regulate gene expression at the transcriptional and post-transcriptional levels in eukaryotes [1]. There are two principal classes of small RNAs in plants, classified according to their biogenesis: microRNAs (miRNAs) and small interfering RNAs (siRNAs) [2, 3]. MicroRNAs are 21–22 nt long and are produced from a single-strand RNA precursor folded into a hairpin. Plant miRNAs are well characterised. They induce post-transcriptional gene silencing (PTGS) principally by triggering messenger RNA (mRNA) degradation, but

they can also induce translational repression [4]. By contrast, siRNAs are 21–24 nt long and are produced from double stranded RNA (dsRNA) precursors resulting from (i) the hybridisation of two complementary RNA strands or (ii) de novo synthesis from a single-stranded RNA as a new complementary strand by RNA-dependent RNA polymerases (RDRs) [5, 6]. siRNAs mediate gene repression at the transcriptional or post-transcriptional level and form several classes differing in terms of predominant size [3].

Transcriptional gene silencing (TGS) takes place in two phases with different specific actors. The pre-establishment phase involves a DNA-dependent RNA polymerase (Pol) II producing aberrant transcripts, RDR6, which produces dsRNAs that are processed by Dicer-like protein 2 (DCL2) and DCL4 to produce 21–22 nt siRNAs that induce DNA methylation by DOMAINS REARRANGED

<figure><img src="image_3.png" /></figure>

BMC

© The Author(s). 2018 **Open Access** This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

<footnote>* Correspondence: stephanie.jaubert@inra.fr
†Clémence Medina and Martine da Rocha are co-first authors and contributed equally to this work.
¹INRA, Université Côte d'Azur, CNRS, ISA, Paris, France
Full list of author information is available at the end of the article</footnote>
METHYLTRANSFERASE (DRM) 2 and, probably, DRM1, through Argonaute 6 - RNA Induced Silencing Complex (AGO6-RISC) and PolIV transcripts [7, 8]. Once RdDM has been established, a stabilisation phase occurs. This phase involves single-strand RNA transcripts produced by PolIV from intergenic or repetitive regions of the genome, which are rendered double-stranded by RDR2 and processed by DCL3 to produce 24 nt hc-siRNAs. Hc-siRNAs are loaded onto AGO4-RISC to initiate RdDM through the DRM2 and, probably, DRM1 proteins at hc-siRNA-generating loci, using other transcripts produced by the DNA-dependent Pol V [8]. One of the main functions of RdDM is maintaining genome integrity, by ensuring that suppressive levels and types of DNA methylation are maintained at transposable elements. RdDM may also modulate the expression of neighbouring protein-coding genes through the spread of DNA methylation [9–12].

siRNAs are also involved in PTGS. They are loaded onto AGO1 or AGO2-RISC, where they generally either induce the cleavage of target transcripts or prevent their translation [3]. The siRNAs involved in PTGS can be classified into several subfamilies according to the origin of the precursor. Natural antisense siRNAs (nat-siRNAs) are generated by the processing of dsRNA precursors derived from endogenous RNAs with complementary sequences, through the action of DCL4 or DCL2, to generate sRNAs of 21 and 22 nt in length, respectively. Phased, secondary, small interfering RNAs (phasiRNAs) are mostly 21 nt siRNAs derived from a RNA converted to dsRNA by RDR6 and processed by DCL4. One well characterised family of phasiRNAs are *Arabidopsis trans*-acting siRNAs (tasiRNAs) [13, 14]. The production of double-stranded phasiRNA precursors is stimulated by one or more upstream miRNAs, such as *TAS3*-derived tasiRNAs in *A. thaliana*, which are generated from non-coding *TAS3* transcripts by miR390 triggers [15, 16].

In plants, siRNAs have been shown to regulate gene expression in various biological processes, including growth, development [15], cell differentiation [17], and plant responses to abiotic and biotic stresses [18–21]. Various PTGS-inducing siRNAs (21–22 nt) have been shown to be related to plant immunity. Examples include nat-siRNAATGB2 from *A. thaliana*, which is induced by *Pseudomonas syringae* pv. *tomato* and plays a positive role in disease resistance by repressing the pentatricopeptide repeats protein-like gene [22]. The nucleotide-binding leucine-rich repeat (NB-LRR) gene family is widely targeted by secondary siRNAs, and phasiRNAs derived from NBS-LRRs play a key role in regulating plant immunity [14, 23]. For example, in Arabidopsis, miR472 and its RDR6-mediated gene silencing help to modulate both PAMP-triggered (PTI) and effector-triggered (ETI) immunity [24]. TGS-mediating

hc-siRNAs have also emerged as major regulators of plant immunity directing DNA methylation and/or histone modification. A role for hc-siRNA-mediated TGS in plant immunity has also come to light, as fungal elicitors induce alterations in the accumulation of certain hc-siRNAs [12]. Moreover, the RdDM machinery has been shown to be involved in plant responses to several pathogens, including *Botrytis cinerea*, *P. syringae*, and *Agrobacterium tumefaciens* [25–27].

Root-knot nematodes (RKN), *Meloidogyne* spp., are highly polyphagous sedentary plant parasites capable of infesting most crop species [28, 29]. After penetrating host roots, RKN larvae migrate toward the vascular cylinder and reprogram gene expression in several vascular root cells, to induce their development into hypertrophied multinucleate giant feeding cells (GCs) [30]. These GCs are metabolically overactive, and serve as the sole source of the nutrients required for RKN development. The growth of the GCs and divisions of the surrounding cells lead to a root deformation known as a knot or gall. The redifferentiation of vascular cells into GCs results from the extensive reprogramming of gene expression in root cells, in response to RKN signals [31]. The expression of genes encoding proteins involved in metabolism, the cytoskeleton, cell cycle, cell rescue, defence, hormones, cellular communication and cellular transport are modified in galls and giant cells from various plant species [30, 32–35]. We are beginning to decipher the genetic pathways modified in infected roots for the formation of GCs, but little is known about the factors regulating this reprogramming of gene expression.

To date, the role of siRNAs in plant-nematode interactions has been little investigated. Two studies provided a first general overview of the small RNA populations expressed in the early feeding sites induced in *A. thaliana* roots by the RKN *M. javanica* [36] and the beet cyst nematode (CN) *Heterodera schachtii* [37]. The first detailed analysis of plant siRNAs expressed in response to *H. schachtii* identified 125 putative *A. thaliana* siRNAs expressed in root feeding sites named syncytia. The methylome of *A. thaliana* syncytia and the associated population of 24 nt siRNAs were recently studied and their abundance was found to be associated with the hypermethylation of transposable elements (TEs) and gene promoters [38]. Moreover, an analysis of the length distribution of Arabidopsis sRNAs in the early developing galls induced by *M. javanica* showed markedly larger numbers of 24 nt reads than in uninfected roots and an overexpression of *miR390* and its secondary *TAS3*-siRNA in galls [36]. Resistance to RKN of Arabidopsis mutant lines for *TAS3a* confirmed the role of these small RNAs in the plant-RKN interaction via the control of Auxin Responsive Factor 3 expression. Overall, these analyses suggest that plant parasitic nematodes use both
microRNAs and the siRNA pathway to manipulate host gene expression at the transcriptional and/or post-transcriptional levels.

We developed a sequencing strategy to characterise the siRNA populations expressed in galls at two key points in giant cell/gall development: 7 days post infection (dpi), corresponding to the phase of successive nuclear divisions without cytokinesis; and 14 dpi, corresponding to the phase of isotropic growth and an increase in DNA levels through endoreduplication without nuclear divisions [39]. In this study, we focused on genomic regions named clusters corresponding to accumulations of siRNAs differentially expressed (DE) between galls and uninfected roots. This analysis provides insight into the loci targeted by siRNAs during the plant-nematode interaction (coding genes, promoting regions, transposable elements). We then identified differentially expressed genes corresponding to differentially expressed clusters with inversely correlated expression profiles, providing biological support for the regulation of these genes by siRNA pathways.

## Results

### Identification of predicted siRNA clusters in galls and uninfected roots

We sequenced small RNAs expressed in uninfected root inter-nodes and galls 7 and 14 dpi induced by the RKN *M. incognita*. Twelve small RNA libraries corresponding to three independent replicates of galls (G) at 7 and 14 dpi and the corresponding uninfected roots (R) were sequenced. The pipeline used for analysis of small RNAs is presented in Fig. 1. Roots galls are composed of nematode and plant tissue. The vast majority of the reads aligned with *A. thaliana* genome (90–95%) with only a small proportion (5–10%) of the reads aligned with *M. incognita* genome [40].

We predicted siRNA clusters, by annotation with ShortStack package [41] for each of the 12 libraries. Most of the sequenced reads mapped on *A. thaliana* genome are located within clusters (Fig. 2a). siRNAs result of the processing by Dicer of a long dsRNA precursor [42]. A DicerCall (DC) score is attributed by ShortStack to each cluster for which the majority of reads are between 20 and 24 nt in size, a characteristic typical of Dicer processing. A DC cutoff of 0.8 (meaning that at least 80% of the reads are between 20 and 24 nt long) was used to distinguish between non-DCL-derived and DCL-derived loci and to exclude all small RNAs corresponding to the degradation products of long RNAs (mRNAs, rRNAs, tRNAs) [43]. Most of the sequenced reads that are located within cluster on *A. thaliana* genome belong to DC-clusters (Fig. 2a, Additional file 1: Figure S1). Moreover, most of the clusters predicted were found to be DC-clusters: about 80%

of the clusters in galls at 7 and 14 dpi, 72% in uninfected roots at 7 dpi and only 58% in uninfected roots at 14 dpi (Fig. 2a, Additional file 1: Figure S1). The lower proportion of DC-clusters in uninfected roots at 14 dpi principally reflects the effect of the R3–14 dpi library, which contained only 44% DC clusters (Additional file 2: Table S1). Only the R3–14 dpi library had a smaller number of reads mapping to DC-clusters (678,080) than the other libraries suggesting that the R3–14 dpi library may contain more degradation products than the other libraries. siRNA clusters produce a mixture of RNAs of different sizes. The DC score (21, 22, 23 or 24 nt) therefore corresponds to the most frequent size of reads mapped to the cluster and does not always reflect the diversity of reads within the cluster (Additional file 1: Figure S1). In all libraries, the majority of DC clusters (between 94 and 97%) had scores of 23 or 24 nt (Fig. 2). Only the DC clusters were retained for further analysis. As previously described [43], we assigned the predicted siRNA DC-clusters into two categories: 20–21–22 nt clusters corresponding to siRNAs and 23–24 nt clusters corresponding to hc-siRNAs.

About twice as many DC clusters were identified in galls as in uninfected roots, with a mean of 67,630 clusters at 7 dpi and 72,556 clusters at 14 dpi for galls, and a mean of 27,151 clusters at 7 dpi and 17,378 clusters at 14 dpi for uninfected roots (Fig. 2a). The mean size of uninfected root clusters was 247 nt at 7 dpi and 202 nt at 14 dpi, whereas the corresponding sizes in galls were 355 nt at 7 dpi and 340 nt at 14 dpi (Fig. 2a; Additional file 2: Table S1). The larger size of clusters in galls indicates that the over-representation of clusters in galls is not due to a prediction bias for size fragmentation. In addition to the much larger number of clusters in galls than in uninfected roots, the coverage of clusters was also clearly different at 7 dpi, with mean coverage rates of 51 and 96 reads per cluster in galls and roots, respectively. At 14 dpi, no strong difference was observed, with a mean coverage of 52 and 60 reads per cluster in galls and roots, respectively.

### Localization of siRNA clusters that are differentially expressed in galls

As siRNA clusters were predicted independently for each sample, the three libraries for each condition (gall or root) did not yield the same lists of predicted siRNA clusters (Additional file 2: Table S1). For the statistical analysis of our data, we constructed a reference set of predicted siRNA clusters with a strategy similar to that used for *Physcomitrella patens* [43]. For each time point, predicted *A. thaliana* siRNA clusters common to at least two of the three libraries in at least one condition (gall or root) were pooled in a reference set of predicted siRNA clusters (Fig. 1). This reference set comprised
<figure><img src="image_4.png" /><figcaption>**Fig. 1** Pipeline of siRNA analysis from prediction to statistical analysis. Data obtained from the SOLiD sequencing of the 12 libraries were cleaned of adaptors and special sequences (snRNA, snoRNA, mitRNA, tRNA, miRNA and pre-miRNA). The Shortstack algorithm mapped and identified clusters corresponding to genomic regions accumulating siRNAs. The algorithm was first run for each library independently. Bedtools was used to identify clusters that were present in at least two out of three libraries in at least one condition (galls or roots). If these clusters were separated by a distance of less than 2 nucleotides, they were then merged (default parameter) and selected to build a reference set of clusters to perform counting and statistical analysis. From counting data, a DESeq statistical analysis was performed to identify clusters that were differentially expressed between gall and root conditions. Only clusters with a DicerCall (DC-clusters) and a minimum coverage of 2 rpmm in all replicates in at least one condition (galls or roots) were considered as biologically relevant</figcaption></figure>

86,264 and 92,578 siRNA clusters from gall and/or root libraries at 7 and 14 dpi, respectively. We compared the expression levels of these predicted siRNA clusters between galls and roots at each time point, by DESeq statistical analysis [44] with an adjusted *P*-value (*p*<sub>adj</sub>) < 0.05 considered significant. Only differentially expressed (DE) clusters supported by a DC score and displaying

significant levels of expression (at least 2 reads per million mapped reads (RPMM)) in the three libraries, for at least one condition (gall or root) were considered to be robust predicted siRNA clusters and were retained for further analysis (Additional file 3: Table S2). An analysis of the genomic positions of predicted differentially expressed siRNA clusters showed that, at 7 dpi, 2871
a

<table><thead><tr><th></th><th>G7</th><th>R7</th><th>G14</th><th>R14</th></tr></thead><tbody><tr><td>1 SOLiD translated reads mapped</td><td>6 551 286</td><td>6 092 737</td><td>6 464 630</td><td>5 119 251</td></tr><tr><td>2 Reads mapped on clusters</td><td>5 982 429</td><td>5 813 141</td><td>5 940 214</td><td>4 821 464</td></tr><tr><td>3 Clusters</td><td>82 343</td><td>37 865</td><td>90 314</td><td>30 325</td></tr><tr><td>4 Reads mapped on DC-clusters</td><td>3 348 679</td><td>2 649 497</td><td>3 926 581</td><td>1 057 681</td></tr><tr><td>5 DC-clusters</td><td>67 630</td><td>27 152</td><td>72 556</td><td>17 378</td></tr><tr><td>6 DC-clusters/clusters (%)</td><td>82</td><td>72</td><td>80</td><td>58</td></tr><tr><td>7 DC-cluster size (nt)</td><td>355</td><td>247</td><td>340</td><td>202</td></tr><tr><td>8 Reads by DC-cluster</td><td>51</td><td>96</td><td>52</td><td>60</td></tr></tbody></table>

b

<figure><img src="image_5.png" /><figcaption>**Fig. 2** Characteristics of the clusters in galls (G) and uninfected roots (R) libraries at 7 and 14 dpi. **a** Average characteristics for the three libraries of each condition (G7, R7, G14 and R14): (1) translated SOLiD color spaced reads mapped on a *A. thaliana* genome; (2) reads mapped on the genomic regions associated with siRNA accumulation (clusters) on *A. thaliana* genome; (3) mean number of clusters per condition; (4) reads mapped on clusters with a Dicer Call (DC-clusters) i.e. when more than 80% of reads on the cluster have a size between 20 and 24 nt; (5) mean number and (6) proportion of DC-clusters; (7) DC-cluster mean size (nt) and (8) mean abundance of reads by DC-cluster. The data for each library are presented in Additional file 2: Table S1. **b** Stacked bar charts of the number of the various categories of DC-clusters (20, 21, 22, 23, 24 nt) for each condition. The proportions of each category are also presented</figcaption></figure>

and 3157 clusters in galls were located within the body of the gene and in the promoter region, respectively (Fig. 3; Additional file 4: Table S3). We found that 2029 and 2295 of the clusters differentially expressed in galls at 14 dpi were located within gene sequence and putative promoter regions, respectively (Fig. 3; Additional file 4: Table S3). The numbers of clusters in putative promoter regions and within the body of genes were therefore similar. Most of these clusters were covered by reads in both gall and root samples. Overall, 275 and 809 clusters upstream from genes at 7 and 14 dpi, respectively, and 328 and 653 clusters within genes at 7 and 14 dpi, respectively, were covered exclusively by reads in a single condition, either galls or roots (Additional file 5: Table S4). Most of the predicted DE clusters were found to be upregulated in galls (Fig. 3; Additional file 4: Table S3). At 7 dpi, the clusters upregulated in galls accounted for 72.0% of the clusters located in promoter regions and 89.7% of those located

within the body of the gene. At 14 dpi, the clusters upregulated in galls accounted for 99.3% of the clusters located in promoter regions and 99.7% of those located within the body of the gene. Only two clusters were among the 20 differentially expressed clusters with the highest fold change in expression, both at 7 and 14 dpi: one of these clusters is located within a gene encoding a lecithin cholesterol acyltransferase (AT3G44830), and the other is located upstream from a gene encoding a homeodomain-like protein (AT2G13960).

The number of clusters was always larger in galls than in uninfected roots (Additional file 6: Table S5) and 23–24 nt was the main DC category of DE clusters regardless of the type of sample (galls or uninfected roots) or the genomic location (in promoter or in gene body) in comparison to the 20–22 nt DC category (Additional file 6: Table S5). However, these two DC categories clearly differed between galls and roots: the number of 23–24 nt clusters in galls was always larger
<figure><img src="image_6.png" /></figure>

<table><thead><tr><th>DE DC-clusters</th><th>7 dpi</th><th>14 dpi</th></tr></thead><tbody><tr><td>In gene body</td><td>2,871</td><td>2,029</td></tr><tr><td>upregulated</td><td>2,574</td><td>2,022</td></tr><tr><td>In promoter region</td><td>3,157</td><td>2,295</td></tr><tr><td>upregulated</td><td>2,274</td><td>2,279</td></tr></tbody></table>

**Fig. 3** Different categories of DC-clusters according to their genomic location. **a** two classes of DC-clusters differentially expressed in galls vs roots were defined: clusters located within gene body or in 2 kb promoter region. **b** Table with the total number of clusters differentially expressed in galls in comparison to uninfected roots at 7 and 14 dpi according to their genomic location. The number of upregulated clusters is indicated

than in uninfected roots, but the opposite pattern was observed for the 20–22 nt category, for which the number of clusters in uninfected roots was always greater than in galls (Fig. 4 and Additional file 6: Table S5). The over-representation differentially expressed clusters in galls appeared therefore to be specific to the 23–24 nt category. This size corresponds to the hc-siRNAs that are known to repress gene expression by targeting transposable elements located in their promoter region and by inducing RdDM. We therefore focused further analyses on the differentially expressed clusters located in putative promoter regions by i) strengthening cluster bioinformatic predictions with transcriptomic data and ii) investigating the presence of transposable elements or repeats derived from transposable elements within these biologically relevant clusters.

## Analysis of differentially expressed clusters of hc-siRNAs and their putative targeted genes in galls

We restricted our analysis to biologically relevant genes, by focusing on predicted DE siRNA clusters located in the putative promoters of genes differentially expressed in galls with inversely correlated expression profiles. The expression of genes located within 2 kb downstream from the siRNA clusters was investigated within DNA chips data listed in NEMATIC [45]. Only a small fraction of the predicted siRNAs targeted genes differentially expressed

in galls (Additional file 3: Table S2). At 7 dpi, 118 genes covered by 129 differentially expressed siRNA clusters and located within the putative promoter region were differentially expressed in galls and the expression profiles of 59 of these genes were inversely correlated with those of the corresponding siRNA clusters, with 37 genes repressed and 22 upregulated in galls (Table 1; Additional file 7: Table S6). Of these 59 DE genes, 14 genes (23.7%) were specifically differentially expressed at 7 dpi and 31 genes (52.5%) were differentially expressed at 7, 14 and 21 dpi (Table 1). Twenty-nine of these 31 genes had similar patterns of expression at these three time points: 20 were downregulated and nine were upregulated in galls, throughout gall development.

At 14 dpi, 116 genes covered by 129 DE siRNA clusters in their promoter regions were differentially expressed in galls, and the expression profiles of 69 of these genes were inversely correlated with those of the siRNA clusters; all of these genes were repressed in galls (Table 2; Additional file 8: Table S7). Fifteen (21.7%) of the 69 DE genes with expression patterns inversely correlated with that of the DE clusters in their putative promoter regions were differentially expressed specifically at 14 dpi, 16 genes (23.2%) were differentially expressed at 7, 14 and 21 dpi. 34 genes (49.3%) were differentially expressed at 14 and 21 dpi, but not at 7 dpi (Additional file 8: Table S7).
<figure><img src="image_7.png" /><figcaption>**Fig. 4** Differentially expressed DC-clusters in galls and roots according to their DC size and their genomic location. Number of differentially expressed (DE) DC-clusters in the twelve galls (G) and roots (R) libraries at 7 and 14 dpi according to their DC size and their genomic location: within gene body or promoter regions. a Number of clusters of the 23–24 nt category. b Number of clusters of the 20–22 nt category</figcaption></figure>

A comparison of the results presented in Additional file 7: Table S6 and Additional file 8: Table S7 identified 13 genes as differentially expressed with expression patterns inversely correlated with those of the DE siRNA-clusters targeting their putative promoters at 7 and 14 dpi (Additional file 9: Table S8). These genes included several with molecular functions relating to catalytic activity (2-oxoglutarate (2OG) and Fe(II)-dependent oxygenase, choline kinase 3, Phosphoglycerate mutase, beta-glucosidase and NADP-dependent malic enzyme 2), or encoding DNA-binding proteins (protein-coding NAC domain-containing protein 58).

## Colocalisation of hc-siRNAs clusters and repeats

Genes targeted by hc-siRNA and RdDM have TEs, or repeats derived from TEs, within their promoter sequences. These sequences are targeted by hc-siRNA, to induce TGS of this region through the induction of RdDM. We investigated the presence of transposon-derived repeats within or in the vicinity of the differentially expressed siRNA clusters. We retrieved the 2 kb upstream from the 5'UTR of DE genes with siRNA clusters in their promoter region and inversely correlated pattern of expression for the gene and the siRNA cluster and investigated the presence of repeats within or close to the cluster by comparing these sequences to Repbase sequences with CENSOR
<table><caption>**Table 1** Top ten most upregulated and downregulated siRNA DC-clusters in galls compared to uninfected roots at 7 dpi</caption><thead><tr><th colspan="3">DC-clusters</th><th colspan="3">Downstream genes</th></tr><tr><th rowspan="2">Cluster name</th><th rowspan="2">variation G/R 7dpi</th><th rowspan="2">gene name</th><th rowspan="2">Gene description (TAIR)</th><th colspan="3">Log2 G/R</th></tr><tr><th>7 dpi</th><th>14 dpi</th><th>21 dpi</th></tr></thead><tbody><tr><td>chr4:5373687-5,374,187</td><td>up</td><td>At4g08455</td><td>BTB/POZ domain-containing protein</td><td>-0.9</td><td>No</td><td>-0.9</td></tr><tr><td>chr1:17450567-17,451,109</td><td>up</td><td>At1g47530</td><td>MATE efflux family protein</td><td>-0.8</td><td>No</td><td>-1.0</td></tr><tr><td>chr3:10387344-10,387,729</td><td>up</td><td>At3g27970</td><td>Exonuclease family protein</td><td>-1.1</td><td>No</td><td>-0.9</td></tr><tr><td>chr2:11040438-11,040,725</td><td>up</td><td>At2g25900</td><td>A. thaliana TANDEM ZINC FINGER PROTEIN 1 (ATCTH)</td><td>-0.8</td><td>-0.8</td><td>-1.4</td></tr><tr><td>chr3:9814548-9,815,346</td><td>up</td><td>At3g26720</td><td>glycosyl hydrolase family 38 protein</td><td>-1.0</td><td>No</td><td>No</td></tr><tr><td>chr3:22340217-22,340,552</td><td>up</td><td>At3g60450</td><td>Phosphoglycerate mutase family protein</td><td>-1.2</td><td>-1.4</td><td>-1.5</td></tr><tr><td>chr3:15356722-15,357,007</td><td>up</td><td>At3g43430</td><td>RING/U-box superfamily protein (zinc finger)</td><td>-0.7</td><td>-0.9</td><td>No</td></tr><tr><td>chr1:25745346-25,745,737</td><td>up</td><td>At1g68570</td><td>membrane localized GA transporter (ATNPF3.1)</td><td>-1.2</td><td>-1.5</td><td>-1.3</td></tr><tr><td>chr1:24996663-24,997,231</td><td>up</td><td>At1g66970</td><td>GLYCEROPHOSPHODIESTER PHOSPHODIESTERASE (GDPD) LIKE 1</td><td>-0.9</td><td>-1.5</td><td>-0.9</td></tr><tr><td>chr2:1681688-1,682,004</td><td>up</td><td>At2g04795</td><td>unknown protein</td><td>-1.2</td><td>No</td><td>No</td></tr><tr><td>chr3:10284167-10,284,467</td><td>down</td><td>At3g27740</td><td>carbamoyl phosphate synthetase (CPS) small subunit (carA)</td><td>0.7</td><td>No</td><td>No</td></tr><tr><td>chr1:26139675-26,140,148</td><td>down</td><td>At1g69530</td><td>Expansin A1</td><td>1.3</td><td>2.9</td><td>3.0</td></tr><tr><td>chr1:10506341-10,507,511</td><td>down</td><td>At1g29980</td><td>choice-of-anchor C domain protein</td><td>1.1</td><td>1.9</td><td>1.4</td></tr><tr><td>chr4:2442915-2,443,871</td><td>down</td><td>At4g04830</td><td>methionine sulfoxide reductase B5</td><td>1.9</td><td>No</td><td>No</td></tr><tr><td>chr3:22214692-22,214,777</td><td>down</td><td>At3g60140</td><td>protein similar to beta-glucosidase</td><td>0.7</td><td>-2.1</td><td>-1.2</td></tr><tr><td>chr2:1281392-1,281,621</td><td>down</td><td>At2g04030</td><td>chloroplast-targeted 90-kDa heat shock protein (CR88)</td><td>1.0</td><td>0.9</td><td>No</td></tr><tr><td>chr1:3322776-3,323,550</td><td>down</td><td>At1g10140</td><td>unknown protein</td><td>0.9</td><td>No</td><td>0.7</td></tr><tr><td>chr5:19447678-19,447,954</td><td>down</td><td>At5g48000</td><td>Member of the CYP708A family of cytochrome P450 enzymes (CYP708A2)</td><td>0.8</td><td>No</td><td>No</td></tr><tr><td>chr2:14572190-14,572,294</td><td>down</td><td>At2g34590</td><td>transketolase family protein</td><td>0.8</td><td>1.1</td><td>0.7</td></tr></tbody></table>

These DC-clusters were all shared by galls and roots and located in promoter regions and with expression patterns inversely correlated with those of the downstream differentially expressed genes. The upregulation (up) or downregulation (down) of the cluster in galls (G) compared to uninfected roots (R) at 7 dpi (variation G/R 7 dpi), the AGI gene name, the description of the encoded protein from TAIR and the log2 values of galls/roots at 7, 14 and 21 dpi obtained from microarrays [32, 45] were indicated

<table><caption>**Table 2** Top ten most upregulated siRNA DC-clusters in galls compared to uninfected roots at 14 dpi</caption><thead><tr><th colspan="3">DC-clusters</th><th colspan="3">Downstream genes</th></tr><tr><th rowspan="2">Cluster name</th><th rowspan="2">variation G/R 14 dpi</th><th rowspan="2">gene name</th><th rowspan="2">Gene description (TAIR)</th><th colspan="3">Log2 G/R</th></tr><tr><th>7 dpi</th><th>14 dpi</th><th>21 dpi</th></tr></thead><tbody><tr><td>chr1:25745518-25,745,718</td><td>up</td><td>At1g68570</td><td>Membrane localized GA transporter (NPF3.1)</td><td>-1.2</td><td>-1.5</td><td>-1.3</td></tr><tr><td>chr1:9412541-9,412,703</td><td>up</td><td>At1g27100</td><td>Actin crosslinking protein</td><td>-1.2</td><td>-0.9</td><td>-1.0</td></tr><tr><td>chr2:18703172-18,703,490</td><td>up</td><td>At2g45380</td><td>Myeloid leukemia factor</td><td>No</td><td>-1.1</td><td>-1.6</td></tr><tr><td>chr4:1063076-1,063,502</td><td>up</td><td>At4g02410</td><td>Concanavalin A-like lectin protein kinase family</td><td>No</td><td>-0.7</td><td>-0.6</td></tr><tr><td>chr5:16944206-16,944,689</td><td>up</td><td>At5g42380</td><td>Calmodulin like 37 (CML37)</td><td>-1.4</td><td>-2.0</td><td>-2.4</td></tr><tr><td>chr3:17722355-17,722,693</td><td>up</td><td>At3g48020</td><td>Hypothetical protein</td><td>No</td><td>-0.8</td><td>-0.8</td></tr><tr><td>chr5:14896137-14,896,289</td><td>up</td><td>At5g37500</td><td>Gated outwardly-rectifying K+ channel</td><td>No</td><td>-0.7</td><td>-0.6</td></tr><tr><td>chr4:8677540-8,679,013</td><td>up</td><td>At4g15230</td><td>Pleiotropic drug resistance 2 (PDR2); ATPase</td><td>No</td><td>-1.0</td><td>-0.7</td></tr><tr><td>chr1:23262107-23,262,634</td><td>up</td><td>At1g62810</td><td>Copper amineoxidase 1 (CuAO1)</td><td>No</td><td>-1.3</td><td>No</td></tr><tr><td>chr4:7855461-7,857,189</td><td>up</td><td>At4g13510</td><td>Plasma membrane localized ammonium transporter (AMT1;1)</td><td>-0.7</td><td>-1.0</td><td>-0.6</td></tr></tbody></table>

These DC-clusters were all shared by galls and roots and located in promoter regions and with expression patterns inversely correlated with those of the downstream differentially expressed genes. The upregulation (up) of the cluster in galls (G) compared to uninfected roots (R) at 14 dpi (variation G/R 14 dpi), the AGI gene name, the description of the encoded protein from TAIR and the log2 values of galls/roots at 7, 14 and 21 dpi obtained from microarrays [32, 45] were indicated
[46]. At 7 dpi, 21 of the 64 clusters differentially expressed in galls, located in promoter regions and with expression patterns inversely correlated with those of the associated DE genes had sequences displaying identity to the *A. thaliana* TE (Fig. 5; Table 3). At 14 dpi, 13 of the 76 clusters differentially expressed in galls, located in promoter regions and displaying expression patterns inversely correlated with those of the associated DE genes displayed had sequences displaying identity to the *A. thaliana* TE (Fig. 5; Table 4). The other DC-clusters displayed no direct sequence identity but were located in the vicinity of repeats. Eight genes differentially expressed at 7 dpi and 17 genes differentially expressed at 14 dpi had siRNA clusters and sequences displaying identity to the sequences of transposable elements in the promoter region, although the transposable elements and siRNA clusters did not overlap (Fig. 5; Table 5; Table 6). These genes correspond to 9 DE clusters at 7 dpi and 20 at 14 dpi that have sequence identity with TE in their vicinity. All together, 46% of the clusters differentially expressed in galls at 7 dpi and 39% of the clusters differentially expressed in galls at 14 dpi display sequence identity with *A. thaliana* TE or are located in vicinity of sequence homologous to TE.

Differentially expressed genes with siRNA clusters colocalised with repeats in their promoter regions overlap with genes differentially expressed in the syncytium induced by cyst nematode with methylation-associated TE profiles

Cyst nematode induce the formation of hypermetallic multinucleate feeding site, named syncytium, that results from the induction of an initial syncytial cell within the root parenchyma that then integrates several hundred of the surrounding cells through cell wall dissolution [47, 48]. Therefore the feeding sites induced by CN and RKN differ by their biogenesis but have similar phenotype and share some molecular pathways [49]. We compared the list of DE genes/DE siRNA clusters in galls displaying inversely correlated patterns of expression and with repeats in the gene promoter to the list of the 526 differentially methylated TE-associated genes corresponding to genes differentially expressed in syncytia. We only identified seven genes as differentially methylated TE-associated genes that were differentially expressed in syncytia [38] and targeted by differentially expressed siRNA clusters in GCs at 7 or 14 dpi (Table 7). These genes were differentially expressed in galls and syncytia, with similar expression patterns, had repeats/TEs in

a Cluster in promoter region & in TE

<figure><img src="image_8.png" /></figure>

Cluster in promoter region & in TE vicinity

<figure><img src="image_9.png" /></figure>

b

<table><thead><tr><th>DE DC-clusters</th><th>7 dpi</th><th>14 dpi</th></tr></thead><tbody><tr><td>In promoter region of inversely correlated DE genes</td><td>64</td><td>76</td></tr><tr><td>In promoter region of inversely correlated DE genes and in TE</td><td>21</td><td>13</td></tr><tr><td>In promoter region of inversely correlated DE genes and in vicinity of TE</td><td>9</td><td>20</td></tr></tbody></table>

**Fig. 5** Location of differentially expressed DC-clusters within promoter regions regarding location of transposable elements. **a** Among the differentially expressed (DE) siRNA clusters located in the promoter of genes differentially expressed in galls with an inversely related expression profile, clusters located within TE (left panel) or in the vicinity of TE (right panel) are the best candidates to regulate gene expression in galls by RdDM pathway. **b** Table with the number of clusters differentially expressed in galls in comparison to uninfected roots according to their location compared to transposable elements: within TE or in the vicinity of TE
<table><caption>**Table 3** siRNA DC-clusters differentially expressed in galls at 7dpi with sequence identity to Arabidopsis transposable elements</caption><thead><tr><th colspan="3">DC-clusters</th><th colspan="3">Downstream genes</th><th>TE</th></tr><tr><th>Cluster name</th><th>variation G/R 7dpi</th><th>Gene name</th><th>Gene description (TAIR)</th><th>Log2 G/R 7 dpi</th><th>Similarity</th><th>TE name</th></tr></thead><tbody><tr><td>chr1:10506341-10,507,511</td><td>down</td><td>At1g29980</td><td>Protein of unknown function</td><td>1.1</td><td>0.8863</td><td>BOMZH2</td></tr><tr><td>chr1:11721631-11,721,864</td><td>up</td><td>At1g32450</td><td>Transmembrane nitrate transporter</td><td>-1.0</td><td>0.9127</td><td>BOMZH2</td></tr><tr><td>chr1:13544313-13,545,488</td><td>down</td><td>At1g36180</td><td>acetyl-CoA carboxylase 2 acetyl-CoA carboxylase 2 (ACC2)</td><td>0.8</td><td>0.8685<br>0.7852</td><td>ATHATN1<br>VANDAL16</td></tr><tr><td>chr1:9412110-9,412,695</td><td>up</td><td>At1g27100</td><td>Actin cross-linking protein</td><td>-1.2</td><td>0.8485</td><td>AtSB2</td></tr><tr><td>chr1:9479725-9,480,634</td><td>up</td><td>At1g27290</td><td>unknown protein</td><td>-1.1</td><td>0.8485<br>1.0000</td><td>ATCOPIA3I<br>ATCOPIA3LTR</td></tr><tr><td>chr2:1281392-1,281,621</td><td>down</td><td>At2g04030</td><td>Chaperone protein htpG family</td><td>1.0</td><td>0.9291</td><td>ATREP13</td></tr><tr><td>chr2:14571706-14,572,100</td><td>down</td><td>At2g34590</td><td>Transketolase family protein</td><td>0.8</td><td>0.9466</td><td>ATREP1</td></tr><tr><td>chr2:14572190-14,572,294</td><td>down</td><td>At2g34590</td><td></td><td></td><td>0.9505</td><td>ATREP1</td></tr><tr><td>chr2:15595033-15,595,856</td><td>down</td><td>At2g37120</td><td>S1FA-like DNA-binding protein</td><td>0.8</td><td>0.9756<br>0.9863</td><td>ATMU7<br>ATMU7</td></tr><tr><td>chr2:8939001-8,939,617</td><td>down</td><td>At2g20750</td><td>expansin B1</td><td>0.7</td><td>0.7231</td><td>BRODYAGA1A</td></tr><tr><td>chr3:10284167-10,284,467</td><td>down</td><td>At3g27740</td><td>carbamoyl phosphate synthetase A</td><td>0.7</td><td>0.7131</td><td>TSCL</td></tr><tr><td>chr3:22214692-22,214,777</td><td>down</td><td>At3g60140</td><td>Glycosyl hydrolase superfamily protein</td><td>0.7</td><td>0.8537</td><td>AtSB2</td></tr><tr><td>chr4:2442915-2,443,871</td><td>down</td><td>At4g04830</td><td>methionine sulfoxide reductase B5</td><td>1.9</td><td>0.9746<br>0.9860</td><td>VANDAL3<br>VANDAL3</td></tr><tr><td>chr4:2444411-2,445,805</td><td>down</td><td>At4g04830</td><td></td><td></td><td></td><td></td></tr><tr><td>chr4:5373687-5,374,187</td><td>up</td><td>At4g08455</td><td>BTB/POZ domain-containing protein</td><td>-0.9</td><td>0.6833</td><td>ATLINE1 6</td></tr><tr><td>chr4:6489475-6,490,728</td><td>up</td><td>At4g10500</td><td>2-oxoglutarate (2OG) and Fe(II)-dependent oxygenase superfamily protein</td><td>-0.8</td><td>0.7176<br>0.8434</td><td>ATHAT8<br>ATHATN4</td></tr><tr><td>chr4:9035216-9,035,789</td><td>down</td><td>At4g15930</td><td>Dynein light chain type 1 family (DIN2)</td><td>1.4</td><td>0.9623</td><td>ATTIRTA1</td></tr><tr><td>chr5:16211264-16,212,381</td><td>down</td><td>At5g40480</td><td>embryo defective 3012</td><td>0.7</td><td>0.7356</td><td>ATREP13</td></tr><tr><td>chr5:16944212-16,945,214</td><td>up</td><td>At5g42380</td><td>calmodulin like 37</td><td>-1.4</td><td>0.7117</td><td>TAG1</td></tr><tr><td>chr5:19447678-19,447,954</td><td>down</td><td>At5g48000</td><td>cytochrome P450, family 708, subfamily A</td><td>0.8</td><td>0.9667</td><td>ATRAN</td></tr><tr><td>chr5:3750834-3,753,083</td><td>up</td><td>At5g11670</td><td>NADP-malic enzyme 2</td><td>-1.4</td><td>0.8554</td><td>ATHATN4</td></tr></tbody></table>

siRNA DC-clusters differentially expressed in galls compared to uninfected roots at 7 dpi located in promoter regions and with expression patterns inversely correlated with those of the associated differentially expressed genes and displaying sequence identity to *A. thaliana* transposable elements (TE). The upregulation (up) or down regulation (down) of the cluster in galls (G) compared to uninfected roots (R) at 7 dpi (variation G/R 7 dpi), the AGI gene name, the description of the encoded protein from TAIR and the log2 values of galls/roots at 7 dpi obtained from microarrays [32, 45] were indicated

their putative promoter regions, predicted siRNA clusters differentially expressed in galls with inversely correlated patterns of expression to the gene and were differentially methylated in syncytia and uninfect ed roots. These seven robust candidates for regulation by RdDM during the plant-CN/RKN interactions encode a receptor protein for CEP1 pep tide (AT5G49660) involved in the maintenance organization of cell files or cell morphology in conductive elements, a xylem nitrate transporter (AT1G32450), a protein involved in oxidative stress (methionine sulfoxide reductase B5, AT4G04830), a 20S proteasome subunit (AT2G05840), a tetraspanin (AT2G23810), a protein similar to beta-glucosidase

(AT3G60140) and a protein of unknown function (embryo defective 3012, AT5G40480).

Discussion

The classification of siRNAs is complex, and their study requires a new analysis at genome level for each biological condition analyzed. For this genome-wide analysis, we used the ShortStack bioinformatics tool [41], which i) was developed for plant genomes, ii) predicts de novo areas of the genome in which small RNAs accumulate, named “clusters” and (iii) carries out statistical analyses of the read counts corresponding to these clusters, comparing the gall to root conditions. This algorithm therefore performs clustering analyses rather
<table><caption>**Table 4** siRNA DC-clusters differentially expressed in galls at 14dpi with sequence identity to Arabidopsis transposable elements</caption><thead><tr><th colspan="3">DC-clusters</th><th colspan="3">Downstream genes</th><th colspan="2">TE</th></tr><tr><th rowspan="2">Cluster name</th><th rowspan="2">variation<br>G/R 14 dpi</th><th rowspan="2">Gene name</th><th rowspan="2">Gene description (TAIR)</th><th colspan="2">Log2 G/R</th><th rowspan="2">Similarity</th><th rowspan="2">TE name</th></tr><tr><th>7 dpi</th><th>14 dpi</th></tr></thead><tbody><tr><td>chr1:17018600–17,021,142</td><td>up</td><td>At1g45015</td><td>MD-2-related lipid recognition<br>domain-containing protein</td><td>No</td><td>-1.6</td><td>0.7642</td><td>TNAT2A</td></tr><tr><td>chr2:14492683–14,493,021</td><td>up</td><td>At2g34350</td><td>Nodulin-like / Major Facilitator Superfamily</td><td>No</td><td>-0.7</td><td>0.8050</td><td>ATTIRX1A</td></tr><tr><td>chr2:18083599–18,083,939</td><td>up</td><td>At2g43590</td><td>Chitinase family</td><td>No</td><td>-1.1</td><td>0.7284</td><td>ATTIRX1B</td></tr><tr><td>chr2:5923051–5,923,711</td><td>up</td><td>At2g14080</td><td>Disease resistance protein (TIR-NBS-LRR class)</td><td>No</td><td>-0.7</td><td>0.8189</td><td>ATLINE1A</td></tr><tr><td>chr2:5924097–5,924,304</td><td>up</td><td></td><td></td><td></td><td></td><td>0.6909</td><td>ATCOPIA95 I</td></tr><tr><td>chr3:10719116–10,720,702</td><td>up</td><td>At3g28600</td><td>P-loop containing nucleoside triphosphate hydrolases</td><td>No</td><td>-1.0</td><td>0.7458</td><td>Sadhu7–2</td></tr><tr><td>chr3:11195491–11,195,822</td><td>up</td><td>At3g29250</td><td>NAD(P)-binding Rossmann-fold superfamily protein</td><td>No</td><td>-0.9</td><td>0.7312</td><td>ATHILA4D LTR</td></tr><tr><td>chr3:22215861–22,216,213</td><td>up</td><td>At3g60140</td><td>Protein similar to beta-glucosidase and is a member<br>of glycoside hydrolase family 1 (DIN2, SRG2)</td><td>0.7</td><td>-2.1</td><td>0.9794</td><td>ATREP13</td></tr><tr><td>chr4:6489725–6,490,158</td><td>up</td><td>At4g10500</td><td>Oxidoreductase, 2OG-Fe(II) oxygenase family protein</td><td>-0.8</td><td>2.3</td><td>0.7176<br>0.7826<br>0.7794<br>0.8434</td><td>ATHATN3<br>ATHAT8<br>ATHATN4<br>ATHATN4</td></tr><tr><td>chr4:6490209–6,490,725</td><td>up</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>chr5:20165643–20,165,928</td><td>up</td><td>At5g49660</td><td>Xylem intermixed with phloem 1</td><td>No</td><td>-1.2</td><td>0.7474</td><td>ATHATN3</td></tr><tr><td>chr5:3750834–3,753,121</td><td>up</td><td>At5g11670</td><td>Malic enzyme (NADP-ME2)</td><td>-1.4</td><td>-1.0</td><td>0.8554</td><td>ATHATN4</td></tr><tr><td>chr5:5716875–5,717,705</td><td>up</td><td>At5g17340</td><td>Putative membrane lipoprotein</td><td>No</td><td>-2.1</td><td>0.7153</td><td>AtSB4</td></tr></tbody></table>

siRNA DC-clusters differentially expressed in galls compared to uninfected roots at 14 dpi located in promoter regions and with expression patterns inversely correlated with those of the associated differentially expressed genes and with sequence identity to *A. thaliana* TE. The upregulation (up) of the cluster in galls (G) compared to uninfected roots (R) at 7 dpi (variation G/R 7 dpi), the AGI gene name, the description of the encoded protein from TAIR and the log2 values of galls/roots at 7 and 14 dpi obtained from microarrays [32, 45] were indicated

than comparing the expression levels of single short sequences. In addition to predicting clusters generating siRNAs, ShortStack also provides information about the small-RNA population of each cluster, including the probability that this cluster results from maturation by a

Dicer protein (DicerCall), size distribution, position on the strand, and the most highly represented sequence at the locus. We identified a large number of siRNA clusters in galls and uninfected roots, with the number of clusters predicted in uninfected root libraries much

<table><caption>**Table 5** siRNA DC-clusters differentially expressed in galls at 7dpi located in the vicinity of Arabidopsis transposable elements</caption><thead><tr><th>DC-cluster name</th><th>Gene name</th><th>Gene description</th><th>TE similarity</th><th>TE name</th></tr></thead><tbody><tr><td>chr1:17450567–17,451,109</td><td>At1g47530</td><td>MATE efflux family protein</td><td>1.3333</td><td>ARNOLDY1</td></tr><tr><td>chr1:26139675–26,140,148 chr1:26140251–26,140,500</td><td>At1g69530</td><td>Expansin A1</td><td>2.3000</td><td>RP1_AT</td></tr><tr><td>chr2:1681688–1,682,004</td><td>At2g04795</td><td>hypothetical protein</td><td>1.6000<br>1.6250</td><td>ATHPOGON3<br>ATHPOGON2</td></tr><tr><td>chr3:3587400–3,587,740</td><td>At3g11410</td><td>Phosphatase 2C</td><td>1.9000<br>1.2333</td><td>ATCOPIA75LTR<br>ATHATN4</td></tr><tr><td>chr3:9814548–9,815,346</td><td>At3g26720</td><td>Glycosyl hydrolase family 38 protein</td><td>2.3333<br>1.7273</td><td>ATREP10B<br>ATHPOGON3</td></tr><tr><td>chr3:20629820–20,629,940</td><td>At3g55610</td><td>delta 1-pyrroline-5-carboxylate synthetase B</td><td>2.8333</td><td>RP1_AT</td></tr><tr><td>chr4:5791975–5,792,114</td><td>At4g09030</td><td>arabinogalactan protein (AGP10).</td><td>1.1818</td><td>BOMZH1</td></tr><tr><td>chr2:2232130–2,232,414</td><td>At2g05840</td><td>20S proteasome subunit PAA2</td><td>1.7619<br>1.8333<br>1.6047<br>1.4730</td><td>SIMPLEGUY1<br>TAG2<br>ATMUNX1<br>SIMPLEGUY1</td></tr></tbody></table>

siRNA DC-clusters differentially expressed in galls at 7 dpi located in promoter regions, with expression patterns inversely correlated with those of the associated differentially expressed genes and located in the vicinity of *A. thaliana* TE. The AGI gene name, the description of the encoded protein from TAIR, the TE similarity and TE name were indicated
<table><caption>**Table 6** siRNA DC-clusters differentially expressed in galls at 14dpi located in the vicinity of Arabidopsis transposable elements</caption><thead><tr><th>DC-cluster name</th><th>Gene name</th><th>Gene description</th><th>TE similarity</th><th>TE name</th></tr></thead><tbody><tr><td>chr1:8383930-8,384,530</td><td>At1g23710</td><td>hypothetical protein</td><td>0.7606</td><td>ATHATN9</td></tr><tr><td>chr1:9412105-9,412,396 chr1:9412541-9,412,703</td><td>At1g27100</td><td>Actin cross-linking protein</td><td>0.8485</td><td>AtSB2</td></tr><tr><td>chr1:11737691-11,737,819</td><td>At1g32460</td><td>hypothetical protein</td><td>0.9656</td><td>ATREP3</td></tr><tr><td>chr1:12890229-12,890,337</td><td>At1g35190</td><td>2-oxoglutarate (2OG) and Fe(II)-dependent oxygenase</td><td>0.7529<br>0.8714</td><td>ATLINE1 5<br>ATHATN2</td></tr><tr><td>chr1:24483137-24,483,294</td><td>At1g65820</td><td>microsomal glutathione s-transferase</td><td>0.9008</td><td>ATHPOGON3</td></tr><tr><td>chr2:743734-743,891</td><td>At2g02680</td><td>Cysteine/Histidine-rich C1 domain family</td><td>0.8144<br>0.8421<br>0.7864</td><td>TNAT1A<br>DT1<br>TNAT1A</td></tr><tr><td>chr2:7123544-7,123,869</td><td>At2g16430</td><td>purple acid phosphatase 10</td><td>0.7983</td><td>ATMUNIX1</td></tr><tr><td>chr2:9979300-9,979,500</td><td>At2g23430</td><td>cyclin-dependent kinase inhibitor protein (KRP1)</td><td>0.8657</td><td>ATTIR16T3A</td></tr><tr><td>chr2:10139020-10,139,268</td><td>At2g23810</td><td>Tetraspanin Tet8</td><td>0.7800</td><td>BRODYAGA1A</td></tr><tr><td>chr2:18703172-18,703,490</td><td>At2g45380</td><td>myeloid leukemia factor</td><td>0.9818</td><td>ATLINEIII</td></tr><tr><td>chr3:9964887-9,966,036</td><td>At3g27020</td><td>metal-nicotianamine transporter YSL6</td><td>0.7614</td><td>ATREP7</td></tr><tr><td>chr3:11031230-11,031,528<br>chr3:11030736-11,031,214</td><td>At3g29034</td><td>transmembrane protein</td><td>0.6516</td><td>ATCOPIN_LTR</td></tr><tr><td>chr4:6144430-6,145,149</td><td>At4g09750</td><td>NAD(P)-binding Rossmann-fold superfamily</td><td>0.6691</td><td>ATLINE2</td></tr><tr><td>chr4:6312052-6,312,340</td><td>At4g10110</td><td>RNA recognition motif (RRM)-containing protein</td><td>0.7590<br>0.7882<br>0.8201</td><td>ATMUN1<br>ATLINE1A<br>ATREP1</td></tr><tr><td>chr5:8671661-8,672,311</td><td>At5g25140</td><td>Cytochrome P450</td><td>0.8462</td><td>ATREP2</td></tr><tr><td>chr5:16944206-16,944,689</td><td>At5g42380</td><td>calmodulin like 37</td><td>0.7117</td><td>TAG1</td></tr><tr><td>chr5:25936470-25,936,892</td><td>At5g64900</td><td>Putative 92-aa protein that is the precursor of AtPep1</td><td>0.8824<br>0.8854</td><td>ATREP6<br>ATREP6</td></tr><tr><td>chr5:25935269-25,935,884</td><td></td><td></td><td></td><td></td></tr></tbody></table>

List siRNA DC-clusters differentially expressed in galls at 14 dpi located in promoter regions, with expression patterns inversely correlated with those of the associated differentially expressed genes and located in the vicinity of *A. thaliana* TE. The AGI gene name, the description of the encoded protein from TAIR, the TE similarity and TE name were indicated

lower than that in galls. We investigated whether this larger number of clusters in the galls was due to the simple fragmentation of clusters into smaller clusters. The larger size of the clusters in galls was not consistent with such a prediction bias. The imbalance for 23–24 nt sequences seems to be characteristic of galls and has already been reported for galls collected at 3 dpi [36]. DESeq statistical analysis identified clusters differentially regulated between galls and uninfected roots. We increased the robustness of our data and the stringency of our analysis by selecting only differentially expressed clusters for which prediction was based on i) reproducible results for the various libraries, ii) a significant level of expression and iii) a high probability of resulting from Dicer cleavage.

The hc-siRNAs associated with TGS are mostly targeted to promoter regions, although some examples of hc-siRNAs targeting the body of the gene have been reported [7, 50]. For the study of hc-siRNAs, we therefore chose to focus on siRNA clusters targeting putative

promoter regions. For the selection of clusters with regulatory action supported by biological evidence, we restricted our study to differentially expressed clusters with a DicerCall of 23–24 nt, covering promoter regions or genes differentially expressed in galls in DNA chip analysis and displaying a pattern of expression inversely correlated with that of the corresponding clusters. However, this strategy was highly restrictive, because only 2.7% of the genes with a cluster in their promoter region displayed transcriptomic regulation in NEMATIC DNA chip data. Similarly, only 5.6% of the DE siRNA clusters target a promoter region of a gene that is differentially expressed in galls at 7 and/or 14 dpi. The NEMATIC data were extracted from transcriptomic analyses performed with microarrays [45]. Microarrays are less sensitive and less exhaustive than new sequencing technologies [51]. New high-throughput sequencing-based analyses of gall transcriptomes are, therefore, required, to obtain a more complete view of the transcriptional reprogramming occurring in galls. The proportions of
<table><caption>**Table 7** siRNA DC-clusters differentially expressed in galls at 7 and/or 14dpi colocalised with differentially methylated regions in cyst nematode feeding sites</caption><thead><tr><th rowspan="2">Gene name</th><th colspan="2">Log2 G/R</th><th rowspan="2">Log2 Sync/R</th><th rowspan="2">Gene description</th><th rowspan="2">Clusters</th><th rowspan="2">CENSOR</th></tr><tr><th>7 dpi</th><th>14 dpi</th><th>5 + 15dpi</th></tr></thead><tbody><tr><td colspan="7"><strong>7 dpi DE DC clusters</strong></td></tr><tr><td>At1g32450</td><td>-1.0</td><td>-1.0</td><td>-3.9</td><td>Transmembrane nitrate transporter</td><td>chr1:11721631-11,721,864</td><td>ATREP4</td></tr><tr><td>At2g05840</td><td>0.8</td><td>No</td><td>0.9</td><td>20S proteasome subunit</td><td>ch r2:2232130-2,232,414</td><td>ATMUNIX1& SIMPLEGUY1 & TAG1</td></tr><tr><td>At3g60140</td><td>0.7</td><td>-2.1</td><td>2.7</td><td>protein similar to beta-glucosidase</td><td>chr3:22214692-22,214,777</td><td>AtSB2</td></tr><tr><td>At4g04830</td><td>1.9</td><td>No</td><td>0.5</td><td>methionine sulfoxide reductase B5</td><td>ch r4:2442915-2,443,871<br>ch r4:2444411-2,445,805</td><td>VANDAL3 VANDAL3</td></tr><tr><td>At5g40480</td><td>0.7</td><td>0.7</td><td>1.3</td><td>Embryo defective 3012</td><td>chr5:16211264-16,212,381</td><td>ATREP13</td></tr><tr><td colspan="7"><strong>14 dpi DE DC clusters</strong></td></tr><tr><td>At2g23810</td><td>-1.0</td><td>-0.9</td><td>-3.3</td><td>Tetraspanin 8</td><td>chr2:10139020-10,139,268</td><td>BRODYAGA1A</td></tr><tr><td>At5g49660</td><td>No</td><td>-1.2</td><td>-1.0</td><td>C-terminally encoded peptide receptor 1</td><td>chr5:20165643-20,165,928</td><td>ATHATN3</td></tr></tbody></table>

siRNA DC-clusters differentially expressed (DE) in galls at 14dpi located in the vicinity of Arabidopsis transposable elements. List of siRNA DC-clusters differentially expressed in galls at 7 dpi and/or 14 dpi located in promoter region, with expression patterns inversely correlated with those of the associated differentially expressed genes and that colocalise with differentially methylated regions in cyst nematode feeding sites. The AGI gene name, the description of the encoded protein from TAIR and the log2 values of galls/roots (G/R) at 7 and 14 dpi obtained from microarrays [32, 45] and the log2 values of syncytia/roots (Sync/R) at 5 and 15 dpi obtained from [38] were indicated

upregulated and downregulated siRNA clusters with inversely correlated expression differed between 7 and 14 dpi. At 7 dpi, 58 clusters (65.5%) were upregulated in galls, whereas, at 14 dpi, 100% of DE clusters were upregulated, suggesting a progressive increase in the number of clusters upregulated in galls. An analysis of the expression profiles of DE genes targeted by siRNA clusters in their putative promoters with an inversely correlated expression profile throughout gall development studied in [32] at 7, 14 and 21 dpi showed that most of the genes targeted by siRNA clusters at 7 dpi were differentially expressed at 7, 14 and 21 dpi, whereas most of the genes targeted by siRNA clusters at 14 dpi are differentially expressed at 14 and 21 dpi, but not at 7 dpi. However, these differences do not suggest any hypotheses concerning the action of siRNA, because these proportions are consistent with the general transcriptomic data obtained over the entire period of gall development [32]. Finally, TEs or sequences displaying identity to repeats were found in the putative promoter region of 40% (at 7 dpi) and 50.7% (at 14 dpi) of the DE genes displaying an inverse correlation of expression with the siRNA clusters in the promoter. These genes are good candidates for regulation by hc-siRNA during gall formation.

hc-siRNAs control transposons and gene expression by inducing DNA methylation. We compared the list of hc-siRNA clusters differentially expressed in galls and supported by biological expression data with the genomic regions that were identified as differentially methylated in feeding sites in *A. thaliana* roots induced by the cyst nematode *H. schachtii* [38]. The role of siRNAs in plant

CN nematode interactions was first highlighted by a significant lower rate of infection with *H. schachtii* in the *A. thaliana* mutants *dcl2-1* [37]. We identified five genes at 7 dpi and two genes at 14 dpi with i) promoter regions displaying some sequence identity to transposable elements ii) differential expression in galls and an inverse correlation of expression with DE siRNA clusters located in the promoter region of the genes concerned and iii) promoter regions differentially methylated in CN feeding sites. These genes are, therefore, robust candidates for regulation by hc-siRNA, through RdDM, during the development of feeding sites induced by plant-parasitic nematodes. Recently, two Arabidopsis mutants *rdr2-rdr6* and *dcl2/dcl3/dcl4* deprived of key factors for RdDM were shown to have lower rates of infection with RKN highlighting the importance of siRNAs and RdDM in Arabidopsis-Meloidogyne interaction [52].

## Conclusions

In this work we provide the first analysis of siRNA clusters expressed in root galls induced by the RKN *M. incognita*. We identified several siRNA clusters that are candidates to modify expression of plant genes in response to RKN infection and these regulations are supported by transcriptomic expression data. These candidates should now be investigated in functional analyses, to confirm i) their regulation by RdDM and ii) their role in plant-nematode interactions. siRNAs will constitute a new field of investigation in studies of the molecular mechanisms underlying plant responses to parasitic nematodes.
# Methods

## Biological material, growth conditions and nematode inoculation

Seeds of *A. thaliana* (ecotype Wassilewskija) were surface-sterilised and sown on Murashige and Skoog (Duchefa) medium agar plates (0.5 x MS, 1% sucrose, 0.8% agar, pH 6.4). Plates were kept at 4 °C for two days, and then transferred to a growth chamber (20 °C with 8 h light and 16 h darkness). *M. incognita* strain “Morelos” was multiplied on tomato plants in a growth chamber (25 °C, 16 h light and 8 h darkness). For in vitro nematode infection, J2 larvae were surface-sterilised with HgCl2 (0.01%) and streptomycin (0.7%) as described elsewhere [53]. We inoculated 25-day-old seedlings grown in vitro individually with 200 sterilised J2s each, resuspended in Phytagel. Seven and 14 dpi, galls were dissected from the infected roots by hand. We also dissected internodes from uninfected roots (without apical and lateral root meristems) with the same age than infected roots for use as a negative control. Samples were immediately frozen in liquid nitrogen and stored at -80 °C. Three independent biological replicates were established for each set of conditions.

## Construction and sequencing of small RNA libraries

Total RNA, including small RNAs (less than 200 nt long), was isolated from galls or uninfected roots at 7 and 14 dpi. Approximately 150 galls or internode fragments from uninfected roots were independently ground into powder in liquid nitrogen, with a mortar. Total RNA was extracted from these samples with the miRNeasy Mini Kit (Qiagen), according to the manufacturer's instructions, with three additional washes in RPE buffer. The quality and integrity of the RNA were assessed with a Bioanalyzer (Agilent). Small RNA libraries were generated by ligation, reverse transcription and amplification (11 cycles) from total RNA (2 µg), with the reagents of the NEBNext Small RNA Library Prep Set for SOLiD. Libraries were then quantified with the Bioanalyzer High Sensitivity DNA Kit (Agilent) and sequenced on a SOLiD 5500 wildfire sequencer (Life Technologies) at the Nice-Sophia Antipolis functional genomics platform (France Génomique, IPMC, Sophia Antipolis, France).

## Bioinformatic siRNA analysis

SOLiD colour spaced reads were translated into sequences for implementation in ShortStack application. For each library, adapters were trimmed and reads matching ribosomal RNA, mitochondrial RNA were removed by performing Blast analyses with the sequences listed in the Rfam database [54]. Reads mapping on sequences corresponding to snRNA, snoRNA, mitRNA, tRNA, miRNA and pre-miRNA were removed in order to keep only siRNAs. Genomic loci accumulating

siRNAs were predicted by using the ShortStack version 3.3 algorithm [41, 55]. ShortStack was run on each library separately with default parameters except: “nohp” mode (no miRNA research) and zero mismatches allowed. ShortStack software version 3.3 mapped trimmed and cleaned reads on a virtual concatenated genome composed by *A. thaliana* genome completed with plastidial and mitochondrial genomes (TAIR10.21) and *M. incognita* genome [56]. From alignments, ShortStack identified a list of coordinates of loci accumulating siRNA. Bedtools [57] with Multi-intersectbed -i option was used to find shared clusters in the different libraries. The output was filtered so that only regions that were present in at least 2 (out of 3) libraries for at least one condition either gall or uninfected root were used to serve as the final reference small RNA locus boundaries. Close clusters with a maximum distance of two nucleotides were merged with bedtools option Merge. ShortStack-count mode under default settings was then used to find relative small RNA abundance on this reference list of each library. Reads mapped on multiple loci were counted on each locus. Counts for siRNA accumulating loci from each replicate were used for differential expression analysis with the R package DESeq. Loci accumulating differentially siRNAs between galls and uninfected roots at a 5% false discovery rate (adjusted P value < 0.05; Benjamini Hochberg adjustment) were retrieved. Localisation of DE clusters was established in gene or in putative promoter gene region defined as the 2 kb upstream coding DNA sequence.

Sequences derived from transposable elements were searched within promoter gene region by using CENSOR algorithm [45] that screens query sequences against a reference collection of repeats. The 2 kb upstream 5'UTR of genes were retrieved and analysed by CENSOR with default parameters.

## Additional files

**Additional file 1: Figure S1.** Percentage of reads of 20, 21, 22, 23, 24 nt within each category of DicerCall clusters in gall (G7, G14) and root (R7, R14) libraries at 7 and 14 dpi. (PPTX 122 kb)

**Additional file 2: Table S1.** Characteristics of mapped reads and clusters predicted by Shortstack 3.3 in each of the three gall libraries (G1-G3) and root libraries (R1-R3) at 7 and 14 dpi. (XLSX 12 kb)

**Additional file 3: Table S2.** Selection of clusters differentially expressed (DE) in galls at 7 and 14dpi. (XLSX 12 kb)

**Additional file 4: Table S3.** Lists of DC-clusters differentially expressed in galls at 7 and 14 dpi and located within a gene body or in a promoter region. (XLSX 3559 kb)

**Additional file 5: Table S4.** Lists of DC-clusters exclusively expressed in galls or roots at 7 and 14 dpi and located within a gene body or in a promoter region. (XLSX 720 kb)

**Additional file 6: Table S5.** Number and proportion of differentially expressed (DE) DC-clusters located in a gene body or in a promoter region presented by DC categories. (XLSX 13 kb)
**Additional file 7: Table S6.** List of DC-clusters differentially expressed (DE) in galls versus uninfected roots at 7 dpi, located within promoter region of genes differentially expressed in galls at 7 dpi with inversely correlated expression profiles. (XLSX 23 kb)

**Additional file 8: Table S7** List of DC-clusters differentially expressed (DE) in galls versus uninfected roots at 14 dpi, located within promoter region of genes differentially expressed in galls at 14 dpi with inversely correlated expression profiles. (XLSX 23 kb)

**Additional file 9: Table S8.** List of DC-clusters differentially expressed (DE) in galls versus uninfected root at 7 and 14 dpi, located within promoter region of genes differentially expressed in galls at 7 and 14 dpi with inversely correlated expression profiles. (XLSX 15 kb)

**Abbreviations**

AGO: Argonaute; DC: DicerCall; DCL: Dicer-like protein; DE: Differentially expressed; dpi: Days post-infection; DRM: Domains Rearranged Methyltransferase; dsRNA: Double stranded RNA; ETI: Effector-triggered immunity; GC: Giant feeding cell; hc-siRNAs: Heterochromatic siRNAs; miRNA: MicroRNA; mRNA: Messenger RNA; nat-siRNA: Natural antisense siRNA; nt: Nucleotide; phasiRNA: Phased, secondary, small interfering RNA; Pol: DNA-dependent RNA polymerase; PTGS: Post-transcriptional gene silencing; PTI: PAMP-triggered immunity; RdDM: RNA directed DNA methylation; RDRs: Dependent RNA polymerases; RISC: RNA Induced Silencing Complex; RKN: Root-knot nematodes; RPMM: Reads per million mapped reads; rRNA: Ribosomal RNA; siRNA: Small interfering RNA; sRNA: Small RNA; tasiRNA: *Trans*-acting siRNA; TE: Transposable element; TGS: Transcriptional gene silencing; tRNA: Transfer RNA

**Acknowledgements**

The authors thank Sébastien Cunnac (IRD, UMR IPE, Montpellier, France) for his help with ShortStack algorithm. The authors thank Christine Lelandais (IPS2, Paris-Saclay, France), Nicolas Bouché (JPB, Versailles, France), Carolina Escobar-Lucas, Javier Cabrera and Virginia Ruiz-Ferrer (University of Toledo, Spain) and Lionel Navarro (ENS, Paris, France) for fruitful discussions. The authors thank the genotoul bioinformatics platform Toulouse Midi-Pyrenees (Bioinfo Genotoul) for providing computing and storage resources.

**Funding**

C.M. was supported by "Santé des Plantes et Environnement" INRA department and Provence Alpes Côte d'Azur fellowship. This work was funded by INRA and by the French Government (National Research Agency, ANR) through the 'Investments for the Future' LabEx SIGNALIFE: program reference #ANR-11-LABX-0028-01 and by the Plant-KBBE program NESTOR (ANR-13-KBBE-0003-06). This Work was supported by France Génomique National infrastructure, funded as part of "Investissement d'avenir" program managed by Agence Nationale pour la Recherche (contrat ANR-10-INBS-09). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Availability of data and materials**

The data set generated during the current study are available in the GEO database (http://www.ncbi.nlm.nih.gov/geo/) with the accession no. GSE100498.

**Authors' contributions**

CM and MDR contributed equally to this work. CM performed the research, bioinformatic and biological data analysis, interpretation; MDR supervised and performed bioinformatic data analysis; NM produced the nematodes; MM and AR collected biological samples; KL performed the preliminary bioinformatic analyses of sequencing results; SJP and BF supervised the experiments and data analyses; SJP, BF and PA supervised the research; CM, SJP and BF wrote the paper. All authors have read and approved the manuscript.

**Ethics approval and consent to participate**

Plant GMO contained use statement No. 564 October 17, 2014.

**Consent for publication**

Not applicable.

**Competing interests**

The authors' declare that they have no competing interests.

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Author details**

¹INRA, Université Côte d'Azur, CNRS, ISA, Paris, France. ²UCA Genomix, Institut de Pharmacologie Moléculaire et Cellulaire, CNRS UMR6097, Sophia Antipolis, Nice, France.

Received: 20 April 2018 Accepted: 21 November 2018

Published online: 18 December 2018

**References**

1. Carthew RW, Sontheimer EJ. Origins and mechanisms of miRNAs and siRNAs. Cell. 2009;136:642–55.

2. Axtell MJ. Classification and comparison of small RNAs from plants. Annu Rev Plant Biol. 2013;64:137–59.

3. Borges F, Martienssen RA. The expanding world of small RNAs in plants. Nat Rev Mol Cell Biol. 2015;16:727–41.

4. Rogers K, Chen X. Biogenesis, turnover, and mode of action of plant microRNAs. Plant Cell. 2013;25:2383–99.

5. Brodersen P, Voinnet O. The diversity of RNA silencing pathways in plants. Trends Genet. 2006;22:268–80.

6. Bologna NG, Voinnet O. The diversity, biogenesis, and activities of endogenous silencing small RNAs in Arabidopsis. Annu Rev Plant Biol. 2014;65:473–503.

7. Matzke M, Kanno T, Daxinger L, Huettel B, Matzke AJ. RNA-mediated chromatin-based silencing in plants. Curr Opin Cell Biol. 2009;21:367–76.

8. Deleris A, Halter T, Navarro L. DNA methylation and demethylation in plant immunity. Annu Rev Phytopathol. 2016;54:579–603.

9. Vaistij FE, Jones L, Baulcombe DC. Spreading of RNA targeting and DNA methylation in RNA silencing requires transcription of the target gene and a putative RNA-dependent RNA polymerase. Plant Cell. 2002;14:857–67.

10. Daxinger L, Kanno T, Bucher E, van der Winden J, Naumann U, Matzke AJM, Matzke M. A stepwise pathway for biogenesis of 24-nt secondary siRNAs and spreading of DNA methylation. EMBO J. 2009;28:48–57.

11. Ahmed I, Sarazin A, Bowler C, Colot V, Quesneville H. Genome-wide evidence for local DNA methylation spreading from small RNA-targeted sequences in Arabidopsis. Nucleic Acids Res. 2011;39:6919–31.

12. Dowen RH, Pelizzola M, Schmitz RJ, Lister R, Dowen JM, Nery JR, Dixon JE, Ecker JR. Widespread dynamic DNA methylation in response to biotic stress. Proc Natl Acad Sci U S A. 2012;109:E2183–91.

13. Vazquez F, Vaucheret H, Rajagopalan R, Lepers C, Gasciolli V, Mallory AC, Hilbert JL, Bartel DP, Crété P. Endogenous trans-acting siRNAs regulate the accumulation of Arabidopsis mRNAs. Mol Cell. 2004;16:69–79.

14. Fei Q, Li P, Teng C, Meyers B. Secondary siRNAs from Medicago NBS-LRRs modulated via miRNA–target interactions and their abundances. Plant J. 2015;83:451–65.

15. Peragine A, Yoshikawa M, Wu G, Albrecht HL, Poethig RS. SGS3 and SGS2/ SDE1/RDR6 are required for juvenile development and the production of trans-acting siRNAs in Arabidopsis. Genes Dev. 2004;18:2368–79.

16. Allen E, Xie Z, Gustafson AM, Carrington JC, Yang M, Padgett RW, Steward R, Chen X, Crete P, Picault N, et al. microRNA-directed phasing during trans- acting siRNA biogenesis in plants. Cell. 2005;121:207–21.

17. Knauer S, Holt AL, Rubio-Somoza I, Tucker EJ, Hinze A, Pisch M, Javelle M, Timmermans MC, Tucker MR, Laux T. A Protodermal miR394 signal defines a region of stem cell competence in the Arabidopsis shoot meristem. Dev Cell. 2013;24:125–32.

18. Vaucheret H. Post-transcriptional small RNA pathways in plants: mechanisms and regulations. Genes Dev. 2006;20:759–71.

19. Padmanabhan C, Zhang X, Jin H. Host small RNAs are big contributors to plant innate immunity. Curr Opin Plant Biol. 2009;12:465–72.

20. Khraiwesh B, Zhu JK, Zhu J. Role of miRNAs and siRNAs in biotic and abiotic stress responses of plants. Biochimica et Biophysica Acta (BBA) - Gene Regulatory Mechanisms. 2012;1819:137–48.

21. Matzke MA, Mosher RA. RNA-directed DNA methylation: an epigenetic pathway of increasing complexity. Nat Rev Genet. 2014;15:394–408.

22. Katiyar-Agarwal S, Morgan R, Dahlbeck D, Borsani O, Villegas A Jr, Zhu JK, Staskawicz BJ, Jin H. A pathogen-inducible endogenous siRNA in plant immunity. Proc Natl Acad Sci U S A. 2006;103:18002–7.

23. Fei Q, Xia R, Meyers BC. Phased, secondary, small interfering RNAs in posttranscriptional regulatory networks. Plant Cell. 2013;25:2400–15.
24. Boccara M, Sarazin A, Thiébeauld O, Jay F, Voinnet O, Navarro L, Colot V. The Arabidopsis miR472-RDR6 silencing pathway modulates PAMP- and effector-triggered immunity through the post-transcriptional control of disease resistance genes. PLoS Pathog. 2014;10(1):e1003883.

25. López A, Ramírez V, García-Andrade J, Flors V, Vera P. The RNA Silencing Enzyme RNA Polymerase V Is Required for Plant Immunity (CS Pikaard, Ed.). PLoS Genet. 2011;7(12):e1002434.

26. Yu A, Lepere G, Jay F, Wang J, Bapaume L, Wang Y, Abraham AL, Penterman J, Fischer RL, Voinnet O, et al. Dynamics and biological relevance of DNA demethylation in Arabidopsis antibacterial defense. Proc Natl Acad Sci U S A. 2013;110:2389–94.

27. Gohlke J, Scholz CJ, Kneitz S, Weber D, Fuchs J, Hedrich R, Deeken R. DNA methylation mediated control of gene expression is critical for development of crown gall tumors. PLoS Genet. 2013;9(2):e1003267.

28. Blok VC, Jones JT, Phillips MS, Trudgill DL. Parasitism genes and host range disparities in biotrophic nematodes: the conundrum of polyphagy versus specialisation. BioEssays. 2008;30:249–59.

29. Abad P, Williamson VM. Plant nematode interaction: a sophisticated dialogue. Adv Bot Res. 2010;53:147–92.

30. Favery B, Quentin M, Jaubert-Possamai S, Abad P. Gall-forming root-knot nematodes hijack key plant cellular functions to induce multinucleate and hypertrophied feeding cells. J Insect Physiol. 2016;84:60–9.

31. Escobar C, Barcala M, Cabrera J, Fenoll C. Overview of root-knot nematodes and giant cells. In: Escobar C, Fenoll C, editors. Advances in botanical research: plant nematode interactions, vol. 73. Oxford: Elsevier; 2015. p. 1–32.

32. Jammes F, Lecomte P, de Almeida-Engler J, Bitton F, Martin-Magniette ML, Renou JP, Abad P, Favery B. Genome-wide expression profiling of the host response to root-knot nematode infection in Arabidopsis. Plant J. 2005;44:447–58.

33. Barcala M, García A, Cabrera J, Casson S, Lindsey K, Favery B, García-Casado G, Solano R, Fenoll C, Escobar C. Early transcriptomic events in microdissected Arabidopsis nematode-induced giant cells. Plant J. 2010;61:698–712.

34. Damiani I, Baldacci-Cresf P, Hopkins J, Andrio E, Balzergue S, Lecomte P, Puppo A, Abad P, Favery B, Hérouart D. Plant genes involved in harbouring symbiotic rhizobia or pathogenic nematodes. New Phytol. 2012;194:511–22.

35. Kyndt T, Denil S, Haegeman A, Trooskens G, Bauters L, Van Criekinge W, De Meyer T, Gheysen G. Transcriptional reprogramming by root knot and migratory nematode infection in rice. New Phytol. 2012;196:887–900.

36. Cabrera J, Barcala M, García A, Rio-Machín A, Medina C, Jaubert-Possamai S, Favery B, Maizel A, Ruiz-Ferrer V, Fenoll C, et al. Differentially expressed small RNAs in Arabidopsis galls formed by *Meloidogyne javanica*: a functional role for miR390 and its TAS3-derived tasiRNAs. New Phytol. 2016;209:1625–40.

37. Hewezi T, Howe P, Maier TR, Baum TJ. Arabidopsis small RNAs and their targets during cyst nematode parasitism. Mol Plant-Microbe Interact. 2008; 21:1622–34.

38. Hewezi T, Lane T, Piya S, Rambani A, Rice JH, Staton M. Cyst nematode parasitism induces dynamic changes in the root epigenome. Plant Physiol. 2017;174:405–20.

39. de Almeida Engler J, Gheysen G. Nematode-induced endoreduplication in plant host cells: why and how? Molecular plant-microbe interactions. Mol Plant-Microbe Interact. 2013;26:17–24.

40. Medina C, da Rocha M, Magliano M, Ratpopoulo A, Revel B, Marteu N, Magnone V, Lebrigand K, Cabrera J, Barcala M, Silva AC, Millar A, Escobar C, Abad P, Favery B, Jaubert-Possamai S. Characterization of microRNAs from Arabidopsis galls highlights a role for miR159 in the plant response to the root-knot nematode Meloidogyne incognita. New Phytol. 2017;216:882–96.

41. Axtell MJ. ShortStack: comprehensive annotation and quantification of small RNA genes. RNA. 2013;19:740–51.

42. Meister G, Tuschl T. Mechanisms of gene silencing by double-stranded RNA. Nature. 2004;431:343–9.

43. Coruh C, Cho SH, Shahid S, Liu Q, Wierzbicki A, Axtell MJ. Comprehensive annotation of Physcomitrella patens small RNA loci reveals that the heterochromatic short interfering RNA pathway is largely conserved in land plants. Plant Cell. 2015;27:2148–62.

44. Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11(10):r106.

45. Cabrera J, Bustos R, Favery B, Fenoll C, Escobar CNEMATIC. A simple and versatile tool for the insilico analysis of plant-nematode interactions. Mol Plant Pathol. 2014;15:627–36.

46. Kohany O, Gentles AJ, Hankus L, Jurka J. Annotation, submission and screening of repetitive elements in Repbase: RepbaseSubmitter and Censor. BMC Bioinformatics. 2006;7:474.

47. Sobczak M, Golinowski W. Structure of cyst nematode feeding sites. In: Cell biology of plant nematode parasitism. Berlin: Springer; 2008. p. 153–87.

48. Kyndt T, Vieira P, Gheysen G, de Almeida-Engler J. Nematode feeding sites: unique organs in plant roots. Planta. 2013;238:807–18.

49. Yamaguchi YL, Suzuki R, Cabrera J, Nakagami S, Sagara T, Ejima C, Sano R, Aoki Y, Olmo R, Kurata T, Obayashi T, Demura T, Ishida T, Escobar C, Root-Knot SS. Cyst nematodes activate Procambium-associated genes in Arabidopsis roots. Front Plant Sci. 2017;13(8):1195.

50. Zhang H, Tao Z, Hong H, Chen Z, Wu C, Li X, Xiao J, Wang S. Transposon-derived small RNA is responsible for modified function of WRKY45 locus. Nature Plants. 2016;2:16016.

51. Hurd PJ, Nelson CJ. Advantages of next-generation sequencing versus the microarray in epigenetic research. Briefings in Functional Genomics. 2009;8:174–83.

52. Ruiz-Ferrer V, Cabrera J, Martinez-Argudo I, Artaza H, Fenoll C, Escobar C. Silenced retrotransposons are major rasiRNAs targets in Arabidopsis galls induced by Meloidogyne javanica. Mol Plant Pathol. 2018;19:2431–45.

53. Caillaud MC, Favery B. In vivo imaging of microtubule Organization in Dividing Giant Cell. Methods Mol Biol. 2016;1370:137–44.

54. Nawrocki EP, Burge SW, Bateman A, Daub J, Eberhardt RY, Eddy SR, Floden EW, Gardner PP, Jones TA, Tate J, et al. Rfam 12.0: updates to the RNA families database. Nucleic Acids Res. 2015;43:D130–7.

55. Shahid S, Axtell MJ. Identification and annotation of small RNA genes using ShortStack. Methods. 2014;67:20–7.

56. Blanc-Mathieu R, Perfus-Barbeoch L, Aury JM, Da Rocha M, Gouzy J, Sallet E, Martin-Jimenez C, Bailly-Bechet M, Castagnone-Sereno P, Flot JF, Kozlowski DK, Cazareth J, Couloux A, Da Silva C, Guy J, Kim-Jo YJ, Rancurel C, Schiex T, Abad P, Wincker P, Danchin EGJ. Hybridization and polyploidy enable genomic plasticity without sex in the most devastating plant-parasitic nematodes. PLoS Genet. 2017;13(6):e1006777.

57. Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26:841–2.

Ready to submit your research? Choose BMC and benefit from:

• fast, convenient online submission

• thorough peer review by experienced researchers in your field

• rapid publication on acceptance

• support for research data, including large and complex data types

• gold Open Access which fosters wider collaboration and increased citations

• maximum visibility for your research: over 100M website views per year

At BMC, research is always in progress.

Learn more biomedcentral.com/submissions

BMC