Hairy roots and calli were obtained from Phaseolus vulgaris L. cv. Negro Jamapa adjusting the procedure described by Estrada-Navarrete et al. (2007). Ten days after infection with A. rhizogenes K599, the primary root was removed. Then plants were grown in hydroponics (B&D solution supplemented with 8 mM KNO₃, pH 6.5) for 11 days in a chamber with 16 h of light and 8 h of dark at 28°C until hairy roots collection. This differs from Estrada-Navarrete et al. (2007) protocol (step 20). For 6 days old calli induction and collection, seeds were germinated as Estrada-Navarrete et al. (2007) protocol (steps 1–7). Then, after seedlings were grown in hydroponics for 2 days with B&D solution, bacteria (steps 9–10) was inoculated in the stem (not in the cotyledonary nodes). Seedlings were grown 6 days in hydroponics inside a chamber to preserve humidity. Common bean non-transgenic roots were obtained from germinated seeds (steps 1–7) grown for 4 days in hydroponics. All samples were frozen in liquid nitrogen and stored at −80°C until RNA extraction.

Total RNA was extracted from hairy roots and calli using the Trizol reagent (Invitrogen). Ten micrograms of each sample were prepared for deep sequencing of small RNA libraries. Sequences ranging from 18 to 30 nt were purified for the construction of the libraries. A small RNA library of hairy roots was prepared following Illumina's Small RNA alternative sample preparation protocol v1.5. The small RNA library of calli was constructed according to the Illumina's TruSeq Small RNA Sample Prep Kit. Libraries were Single Read-sequenced using the Genome Analyzer IIx (GAIIx) and the Illumina Cluster Station at the Instituto de Biotecnología (Universidad Nacional Autónoma de México). Degradome library construction for the hairy roots sample was performed as described by Ma et al. (2010) with some changes. Approximately 150 ng of poly(A)⁺ RNA was used to anneal with biotinylated random primers. Strapavidin capture of RNA fragments was performed using biotinylated random primers. A 5′ adaptor ligation was only performed to those RNAs containing 5′-monophosphates. The library was sequenced with the Illumina's Cluster Station and the GAIIx using the 5′ adapter only. This resulted in the sequencing of the first 36 nt of the inserts that represent the 5′ ends of the original RNAs (LC Sciences). All sequence data is available at the European Nucleotide Archive (EMBL-EBI) with the accession number PRJEB7993.

Sequence length distribution and alignment (PatMaN; Prüfer et al., 2008) of small RNAs were performed using the UEA small RNA analysis toolkit (Version 2.5.0; Stocks et al., 2012). Small RNAs (18–26 nt) were aligned to the A. rhizogenes plasmid Ri T-DNA region without mismatches (GenBank accession number: EF433766). Whole T-DNA alignments were seen with the Integrative Genomics Viewer (IGV; Thorvaldsdóttir et al., 2013). RNA secondary structures of the T-DNA loci were predicted with the Vienna RNA secondary structure server using the minimum free energy algorithm (MFE; Hofacker, 2003). P. vulgaris genes targeted by T-DNA derived sRNAs were predicted in silico using the plant small RNA target analysis server psRNATarget (Dai and Zhao, 2011) with a maximum expectation threshold of 0–2. P. vulgaris (Phytozome v9.0) spliced mRNA transcripts without alternative splice variants were used for target prediction analysis. Multiple alignment of rolA sequences from different Ri plasmids (pRi1724/gb: AP002086, pRi8196/gb:M60490, pRiA4/gb:X12579) was performed using ClustalW2 with default settings (Larkin et al., 2007) and viewed with Jalview (Waterhouse et al., 2009). Previously reported microRNA sequences from P. vulgaris were used for the identification of microRNAs in the small RNA libraries and for the non-transgenic root miRNA accumulation analysis (Peláez et al., 2012). Normalized read frequencies (RPM) of miRNAs from the same family were added for miRNA accumulation analysis.

Two strategies were used for the identification of loci producing phasiRNAs of 21 nt in length (Chen et al., 2007; Axtell, 2010). The algorithm developed by Chen et al. (2007) was used through the UEA small RNA analysis toolkit implementation (Figures 5B,C). Only sRNAs derived from the T-DNA and the plasmid Ri T-DNA sequence (gb: EF433766) were used as input. Moreover, this algorithm was performed for the identification of phasing-generating loci from P. vulgaris used as controls (Supplementary Table 9). Combined abundances of two and five bins were used in the methodology developed by Axtell (2010) to calculate the phase ratio of rolB and CUS loci, respectively.

Total RNA was extracted using Trizol reagent (Invitrogen). RNA (20 ug) was separated by 15% PAGE/8 M urea/ 1x TBE buffer. Gels were electro-blotted to a Hybond-N⁺ membrane (GE Healthcare) and UV cross-linked twice. Oligonucleotide probes (Supplementary Table 10) were labeled using [γ-³²P]-ATP (Perkin-Elmer) and T4 PNK (New England Biolabs). Labeled probes were purified with the Quick spin-oligo columns (Roche) before hybridization. Hybridizations were performed at 42°C overnight in UltraHyb-oligo solution (Ambion). After hybridization, membranes were washed twice in 2 × SSC/0.1% SDS and exposed to the Phosphor Screen System (GE Healthcare). The screen was scanned in a Storm 860 Gel and Blot Imaging System (GE Healthcare). As loading control, an oligonucleotide probe complementary to the U6 small nuclear RNA was used. Signal intensities were quantified using the ImageQuant 5.2 software (Molecular Dynamics). The U6 signal intensity in each blot was used for normalization and calculation of expression ratios. The same loading control was employed for membranes used more than once.

Perfect alignments of degradome sequences to the T-DNA were performed using PatMaN (Prüfer et al., 2008). The PAREsnip tool was used for the discovery of small RNA-guided cleavage sites with default parameters (Folkes et al., 2012). The small RNAs derived from the T-DNA, the T-DNA sequence (gb: EF433766) and the degradome sequences were used as input for the characterization of sRNA-mediated cleavages on the A. rhizogenes plasmid Ri T-DNA region. Primary transcripts (Phytozome v9.0) were used as input for the discovery of P. vulgaris targets cleaved by T-DNA derived sRNAs. Degradome sequences from seedlings of common bean were explored in the same way as input and used as control (Formey et al., 2015). Host targets were discarded if other sRNA from the library was also predicted to cleave the transcript and presented a better hit in at least one of the parameters evaluated such as alignment score, category, reads abundance, or p-value. Arabidopsis homolog targets were functionally grouped according to the Gene Ontology (GO) categories (Berardini et al., 2004). GO enrichment analysis was performed through the Plant GeneSet Enrichment Analysis Tool kit under default parameters (Yi et al., 2013). MicroRNAs (Supplementary Table 4) were used for analysis of validated target-sRNA interactions in the degradome library (Supplementary Table 5).

Detection and quantification of transcripts accumulation in hairy roots was performed using the iQ5 Real-Time PCR Detection System and the iQ5 Optical System Software (Bio-Rad). First strand cDNA was synthesized from DNA-free RNA with the RevertAid H Minus First Strand cDNA Synthesis kit (Fermentas). The qRT-PCR reactions were carried out in triplicate for three biological replicates with the Maxima SYBR green-fluorescein qPCR master mix (Fermentas). The gene-specific oligonucleotides belonging to the T-DNA genes were tested also on non-transformed tissues as controls (Supplementary Table 11). The melting temperature was set at 60°C. Quantification was based on a cycle threshold (Ct) value and the genes expression levels were normalized with the elongation factor 1 (EF1) reference gene.

In order to identify small RNAs derived from the T-DNA of A. rhizogenes in the hairy root disease, small RNA libraries from hairy roots, and calli (formed before the emergence of hairy roots at the wound site) were generated. The presence and expression of several genes from the T-DNA in hairy roots (Supplementary Figure 1) were confirmed. Small RNA sequencing yielded 12,396,242 and 27,272,479 total raw reads for the hairy roots and calli samples, respectively. Small RNA sequences (18–26 nt) in hairy roots showed that the 21-nt (11%) and the 24-nt (32%) classes were the most abundant (redundant sequences). The 21-nt class was the second most diverse (unique sequences) in hairy roots and the fourth in A. rhizogenes-induced calli (Figure 1A). Perfect sequence alignment of small RNAs against the T-DNA of A. rhizogenes revealed the presence of several distinct and abundant ArT-sRNAs (18–26 nt) in hairy roots (Figure 1B). In hairy roots, 3176 distinct ArT-sRNAs corresponding to 17,000 small RNA reads per million (RPM) were detected. In contrast, only 563 ArT-sRNAs (unique sequences) corresponding to an abundance of 90 RPM were derived from T-DNA in calli (Figure 1C; Supplementary Table 1; De Paoli et al., 2009). Differences in the amount of ArT-sRNAs detected between hairy roots and calli may be a consequence of the amount of transformed cells in each sample. Most of the ArT-sRNAs observed in hairy roots aligned with transcriptional units of the T-DNA belonging to genes, such as ORF2, ORF8, rolA, rolB, ORF13, ORF14, and the opine synthase gene CUS (Figure 1B; Table 1), with the oncogene rolA being the major source of ArT-sRNAs in hairy roots. Interestingly, virtually no ArT-sRNAs derived from the oncogene rolC were detected in hairy roots. Also, nothing but 134 unique ArT-sRNAs were shared between hairy roots and calli.

Moreover, the impact of ArT-sRNAs production on the accumulation of microRNAs in hairy roots and calli was measured. Members of conserved miRNA families, such as miR156, miR157, miR159, miR160, miR162, miR164, miR166, miR167, miR168, miR169, miR319, miR390, miR393, miR394, miR396, and miR408 were considered. Production of miRNAs in hairy roots was similar compared with non-transgenic roots according to northern blot analysis for eight miRNAs (Figure 2A). Relative miRNA expression estimation using the frequencies of miRNAs detected in the hairy roots and the calli libraries showed that the majority of miRNAs analyzed accumulated less in calli than in hairy roots. Substantial differences in accumulation were detected in these two samples for microRNAs miR394, miR408, and miR319 (Figure 2B). These differences in accumulation were also observed for these three families of microRNAs between the read frequencies of calli and a previously reported small RNA library from non-transgenic roots (Supplementary Figure 2; Peláez et al., 2012). MicroRNAs miR408 and miR394 were up-regulated in hairy roots and non-transgenic roots compared with calli, whereas miR319 was down-regulated in both type of roots. Also, miR393, a microRNA involved in plant defenses and auxin signaling, was more abundant in calli than in hairy roots (Figure 2B; Peláez and Sanchez, 2013). Interestingly, miR319 was the most abundant miRNA detected in calli. This was not the case for small RNA libraries of common bean from roots, leaves, flowers, and seedlings (Peláez et al., 2012).

To unravel the probable biogenesis and action of common bean ArT-sRNAs, their size, strand origin, and phased patterns were analyzed. It was found that ArT-sRNAs were largely small RNAs of 21-nt in length in hairy roots and calli (Figure 3A). In hairy roots, 52.2% of the ArT-sRNAs (unique sequences), representing 67.9% of the total sequences that aligned to the T-DNA, were of 21-nt in length. Only 5% of the total sequences corresponded to the 24-nt class that is involved in transcriptional regulation. The 22-nt class of ArT-sRNAs was the second most diverse (18%) and abundant (19.7%). A similar pattern of accumulation was observed for the ArT-sRNAs detected in calli, which indicates that silencing of T-DNA transcripts may occur at the post-transcriptional level. Also, ArT-sRNAs were produced from the sense and antisense strands, suggesting that ArT-sRNAs may be generated from perfect long double stranded RNAs (Figure 3B). In hairy roots, 57% of redundant sequences of the ArT-sRNAs found aligned to the sense strand while 43% aligned to the antisense strand. Taking a closer look to particular genes encoded in the sense or antisense strand, the proportion of ArT-sRNAs that aligned to either strand was variable (Figure 4). Furthermore, phased secondary 21-nt small interfering RNAs characteristic of PHAS loci were identified for ORF8, rolA, rolB, ORF13, ORF14 and CUS genes (Figures 5A–C). Genes like rolA, ORF8, ORF13, and ORF14 exhibited one sRNA cluster and one predominant phase, whereas rolB and CUS presented two clusters without one main phase. Some phased small RNAs of the rolA transcript were among the most abundant sRNAs detected (Figure 5A).

Forty percent of the ArT-sRNAs of 21-nt in length that were detected in this analysis had a 5′ terminal uracil (Supplementary Figure 3). Small RNAs with this characteristic are recruited by AGO1. To gain insight into the cleavage of T-DNA transcripts through AGOs/ArT-sRNAs associations, a degradome (also called parallel analysis of RNA ends, PARE) library from hairy roots was deep-sequenced and analyzed. Sequence alignment resulted in the identification of 446 distinct sequences that aligned perfectly to the T-DNA of A. rhizogenes. Degradome sequences corresponding to all of the coding sequences of the T-DNA genes except for ORF13 were detected (Supplementary Figures 4, 5). Most of the reads detected had the sequence of the rolA and CUS genes. Discovery analysis of sRNA/target interactions using degradome sequences and ArT-sRNAs from hairy roots exposed 180 interactions between ArT-sRNAs and the transcripts of the T-DNA with significant p-values at 73 unique cleavage positions (Supplementary Table 2). According to degradome read abundances (Addo-Quaye et al., 2009), 5 interactions fell into category 0 (highly likely to occur), 116 into category 2, and 59 into category 4. From almost all of the interactions found, around 95.5%, pointed out to cleavage sites of the rolA transcript, mainly from the 3′ untranslated region (Figure 6A). These results strongly suggest that AGOs/ArT-sRNAs associations occur. Furthermore, one of the ArT-sRNAs that fell into category 0, and thus cleaved the most abundant fragment in the transcript of rolA, was indeed the most abundant ArT-sRNA in the library. This ArT-sRNA of 22-nt in length with 3304 absolute read counts in the library presented another cleavage position including two mismatches (Figure 6A). This cleavage position in the transcript of rolA was consistent with the phasing register of the detected rolA phase (Figure 5A). In addition, the sequence of this ArT-sRNA was found to be highly conserved among rolA sequences of different Ri plasmids (pRi8196, pRi1724, pRiA4, and pRi2659), which suggests functional roles for this small RNA (Figure 6B).

ArT-sRNAs may present the required complementarity base pairing with the transcripts of P. vulgaris to cleave them. To search for complementarity base pairing between ArT-sRNAs and host transcripts (P. vulgaris), a target prediction analysis was performed using the plant small RNA target analysis server psRNATarget (Dai and Zhao, 2011). The identification of 399 unique ArT-sRNAs/target interactions using a stringent cut-off threshold of the complementarity score (Supplementary Table 3) was determined. The interactions detected were composed by 299 distinct transcripts and 326 distinct ArT-sRNAs. To provide further evidence of ArT-sRNAs-mediated cleavage of host gene transcripts, the degradome library obtained from hairy roots was analyzed in pursuit of putative interactions. MicroRNAs (Supplementary Table 4) were used for analysis of conserved validated target-sRNA interactions in the degradome library as controls. Most of the microRNA-target interactions found corresponded to conserved validated target-microRNA interactions (Supplementary Table 5). In total, 402 interactions between host targets and ArT-sRNAs were detected using PAREsnip with default parameters (Supplementary Table 6). The total number of unique host transcripts predicted to be cleaved by 346 unique ArT-sRNAs was 228. No interactions were found between ArT-sRNAs and these 228 transcripts using the degradome sequencing of non-transformed common bean seedlings (Formey et al., 2015). Most of the ArT-sRNAs that define these interactions were of 21-nt in length (53.4%). Taking into account all the interactions, 58 fell into category 0, 50 into category 1, 155 into category 2, 3 into category 3, and 136 into category 4. The host transcript targets identified correspond to a diverse set of protein families (Supplementary Table 7). Gene Ontology (GO) enrichment analysis was performed for the 98 of the unique host targets (Supplementary Table 8). The two most enriched gene sets within the cellular component category were cellular (GO: 0044464, p-value = 1.91E-09, FDR = 5.97E-07) and intracellular components (GO: 0044424, p-value = 2.97E-07, FDR = 6.96-05). Metabolic (GO: 0008152, p-value = 8.45E-12, FDR = 7.58E-08) and cellular processes (GO: 0009987, p-value = 7.68E-11, FDR = 3.44E-07) were the most enriched gene sets observed from the biological process category. The response to stress category was also detected by the GO enrichment analysis (GO: 0006950, p-value = 2.49E-08, FDR = 4.31E-05) suggesting probable regulation of host transcripts by ArT-sRNAs as a counter-defensive strategy. Given that most of the interactions found by PAREsnip were validated interactions (Folkes et al., 2012) and that we found several interactions between ArT-sRNAs and P. vulgaris transcripts, our findings suggest that ArT-sRNAs could target common bean transcripts.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.