ARTICLE https://doi.org/10.1038/s41467-020-19129-6 OPEN Identification of key sequence features required for microRNA biogenesis in plants Arantxa M. L. Rojas 1, Salvador I. Drusin 1,2, Uciel Chorostecki 1,7,8, Julieta L. Mateos 1,9, Belén Moro 1,10, Nicolas G. Bologna 1,10, Edgardo G. Bresso 1, Arnaldo Schapire1, Rodolfo M. Rasia 1,3, Diego M. Moreno 4,5 & Javier F. Palatnik 1,6 ✉ 1234567890():,; MicroRNAs (miRNAs) are endogenous small RNAs of ∼21 nt that regulate multiple biological pathways in multicellular organisms. They derive from longer transcripts that harbor an imperfect stem-loop structure. In plants, the ribonuclease type III DICER-LIKE1 assisted by accessory proteins cleaves the precursor to release the mature miRNA. Numerous studies highlight the role of the precursor secondary structure during plant miRNA biogenesis; however, little is known about the relevance of the precursor sequence. Here, we analyzed the sequence composition of plant miRNA primary transcripts and found specifically located sequence biases. We show that changes in the identity of specific nucleotides can increase or abolish miRNA biogenesis. Most conspicuously, our analysis revealed that the identity of the nucleotides at unpaired positions of the precursor plays a crucial role during miRNA bio- genesis in Arabidopsis. 1 IBR (Instituto de Biología Molecular y Celular de Rosario), CONICET and Universidad Nacional de Rosario, Rosario 2000, Argentina. 2 Área Física, Departamento de Química-Física, Facultad de Ciencias Bioquímicas y Farmacéuticas, Universidad Nacional de Rosario, S2002LRK Rosario, Santa Fe, Argentina. 3 Área Biofísica, Facultad de Ciencias Bioquímicas y Farmacéuticas, Universidad Nacional de Rosario, Suipacha 531, S2002LRK Rosario, Santa Fe, Argentina. 4 Instituto de Química de Rosario (CONICET-UNR), Suipacha 570, S2002LRK Rosario, Santa Fe, Argentina. 5 Área Química General e Inorgánica, Departamento de Química-Física, Facultad de Ciencias Bioquímicas y Farmacéuticas, Universidad Nacional de Rosario, Suipacha 531, S2002LRK Rosario, Santa Fe, Argentina. 6 Centro de Estudios Interdisciplinarios, Universidad Nacional de Rosario, Rosario 2000, Argentina. 7Present address: Barcelona Supercomputing Centre (BSC-CNS), Barcelona (08034), Spain. 8Present address: Institute for Research in Biomedicine (IRB), The Barcelona Institute of Science and Technology, Barcelona (08028), Spain. 9Present address: Instituto de Fisiología, Biología Molecular y Neurociencias (IFIBYNE), CONICET-UBA, Buenos Aires (1428), Argentina. 10Present address: Centre for Research in Agricultural Genomics (CRAG), CSIC-IRTA-UAB-UB, Campus UAB, Barcelona (08193), Spain. ✉email: palatnik@ibr-conicet.gov.ar NATURE COMMUNICATIONS | (2020)11:5320 | https://doi.org/10.1038/s41467-020-19129-6 | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-19129-6 M icroRNAs (miRNAs) are one class of small RNAs that all the processing steps necessary to release the mature miRNA perform essential functions as post-transcriptional from the precursor are carried out in the nucleus by a complex regulators of gene expression in animals and plants. formed by DICER-LIKE1 (DCL1) and accessory proteins, such as Plant miRNAs are transcribed as longer transcripts by RNA HYPONASTIC LEAVES1 (HYL1) and SERRATE (SE) (reviewed polymerase II, which are capped and polyadenylated (reviewed in in refs. 1,2). DCL1 performs cuts in the foldback at positions refs. 1,2). These primary transcripts harbor an imperfect foldback flanking the miRNA, releasing two paired RNA strands of ~21 nt structure with the miRNA embedded in one of its arms. In plants, with two nt 3′ overhangs, the miRNA/miRNA* duplex. One of these small RNAs, generally the miRNA, becomes incorporated into an ARGONAUTE (AGO) protein where it provides a 5p sequence specificity to identify target RNAs (reviewed in refs. 3,4). Plant miRNA The Arabidopsis genome encodes ten AGO proteins and the precursor 3p miRNA pathway usually recruits AGO1 that preferentially binds small RNAs with 5′ U, which is the most common nucleotide at Nucleotide the 5′-end of plant miRNAs3,4. pairs –18 1 21 38 The precise excision of the mature miRNA from the precursor position Lower miRNA/miRNA* Upper stem duplex stem is a key aspect of miRNA biogenesis, and relevant for the specific regulation of target RNAs. Many evolutionarily conserved miR- b c NAs regulate transcription factors with important roles in plant Arabidopsis thaliana Eudicots development and hormone signaling, hence, mutations in the WC miRNA biogenesis machinery lead to important defects in plant Wobbles MM growth and/or lethality1. Plant miRNA precursors are largely variable in size and shape, yet they are accurately processed to release rather specific miRNA sequences5–7. Previous work has 4 3 identified structural determinants present in the precursors that 30 3 30 guide the miRNA processing machinery to produce specific cuts Percentage Percentage Percentage Percentage 20 20 2 along the precursor sequence5–12. 2 Numerous precursors have a 15–17 bp dsRNA segment above 1 10 1 10 an inner loop and below the miRNA/miRNA* duplex, which 0 0 0 0 guides a first cleavage by DCL1 at the proximal part of the pre- AU/UA GC/CG GU/UG AU/UA GC/CG GU/UG cursor (base-to-loop processing)7,9–11. Other precursors have a CA/AC CU/UC AG/GA UU AA CC GG CA/AC CU/UC AG/GA UU AA CC GG dsRNA segment above the miRNA/miRNA* duplex and below a small terminal loop5,12–14. In these precursors, the processing d e Arabidopsis thaliana Eudicots machinery produces a first cut that releases the terminal loop and First Second First Second then continues towards the base of the precursor to release the Cleavage sites miRNA/miRNA* duplex (loop-to-base processing)5,13,14. In either case, after a first cut is generated, DCL1 produces a sub- Base pair –1 1 2 3 21 22 23 24 position –1 1 2 3 21 22 23 24 sequent cut ~21 bp away from the first one, releasing a miRNA/ 59 90 79 95 77 80 90 80 Paired 57 92 78 86 79 87 87 70 miRNA* duplex7,15. 41 10 21 5 23 20 10 20 Unpaired 43 8 22 14 21 13 11 26 Although previous work has shown the essential role of the 33 18 49 44 57 26 31 59 40 AU/UA 37 36 20 51 39 37 36 36 41 precursor secondary structure for miRNA biogenesis, little is 23 59 30 49 20 46 46 20 36 GC/CG 35 14 55 22 45 42 40 43 22 known about the importance of the precursor primary sequence. 3 13 0 2 0 8 13 2 5 GU/UG 7 7 17 5 2 0 11 8 6 Here, by analyzing precursors from different plant species we 25 7 15 2 2 3 2 8 8 AC/CA 6 24 2 9 4 3 0 2 6 found sequence biases at specific locations. Based on this infor- 2 2 3 2 8 3 2 5 3 CU/UC 4 3 2 7 4 5 3 1 4 mation, we were able to increase or decrease the processing 2 0 2 0 2 2 0 3 1 AG/GA 1 0 1 2 0 2 1 0 1 efficiency by introducing changes in the precursor sequence. 2 0 2 0 8 0 5 2 2 AA 2 2 0 1 0 5 0 2 2 Furthermore, our analysis revealed that the identity of nucleotides 11 2 0 2 3 3 2 2 3 UU 3 11 1 1 4 3 1 0 2 at mismatched positions along the precursor structure can play a 0 0 0 0 0 8 0 0 1 GG 1 0 1 0 0 0 4 0 0 significant role during miRNA processing. Most conspicuously, 0 0 0 0 0 0 0 0 0 CC 0 1 0 1 0 1 0 0 1 we found that C-C mismatches consistently impair miRNA bio- 0 0 0 0 0 0 0 0 0 ssRNA 3 2 1 1 1 2 3 6 9 genesis, and that they are absent from DCL1 cleavage sites. Average Average Percentage Results 0 100 Identification of sequence biases in plant miRNA precursors. Fig. 1 Frequency of nucleotide pairs in plant miRNA precursors. To examine whether MIRNAs from Arabidopsis and other a Schematic representation of a plant microRNA precursor. The microRNA/ angiosperms show any nucleotide bias beyond the miRNA microRNA* is indicated in blue, and the DCL1 cleavage site with dashed sequence, we analyzed miRNA precursors from 30 eudicot species green lines. b, c Overall nucleotide pairs composition in conserved MIRNA (Source Data). First, we determined the frequency of precursors of Arabidopsis thaliana (b) and eudicots (c). Pie chart indicating Watson–Crick base pairs (WC), G-U wobbles, and mismatches the percentage of Watson–Crick pairs (WC), G-U wobbles pairs and (MM) in the secondary structure of Arabidopsis precursors, mismatches (MM). Bars indicate the percentage of each pair. d, e Base pair considering 56 pairs in length, therefore including the miRNA/ composition at DCL1 cleavage sites of miRNA precursors from Arabidopsis miRNA* as well as regions below and above this duplex (Fig. 1a). thaliana (d) and eudicots (e). The first and second cleavage sites are We found that ~65% of the positions were paired, which is indicated. Positions selected for each group of precursors are indicated in consistent with the structured nature of the miRNA precursors Supplementary Fig. 1. Numbers in each square indicate the percentage of (Fig. 1b, c). Still, the G-U wobble (~8%) was less frequent than a each variant on each position. A-U or G-C canonical Watson–Crick base pair (>30% each one) 2 NATURE COMMUNICATIONS | (2020)11:5320 | https://doi.org/10.1038/s41467-020-19129-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-19129-6 ARTICLE (Fig. 1b, c). We found similar proportions when we analyzed the miR172/pri-miR172 was consistently higher for the MIR172 C-G MIRNAs in different eudicots (Fig. 1b, c). mutant variant compared with the wt MIR172 C-U (Fig. 2f, i). Despite mismatches being less frequent than paired positions, The MIR172 C-G variant accumulated slightly more miR172* we were able to analyze 10273 mismatched positions belonging to than the wt precursor, but the difference was not statistically 1681 precursors of 30 eudicot species. Interestingly, we observed significant (Supplementary Fig. 4b). variation in the identity of the nucleotides composing the The change C-U to C-G is located in a structured region of the mismatches. The most common pair was U-U present in ~3% MIR172A precursor (Fig. 2a). To test whether this single change of the positions analyzed. On the other hand, C-C and G-G affects the structure of the region, we probed the precursor using mismatches were present in <0.5% of the cases (Fig. 1b, c). Terbium (Tb III), which induces the cleavage of RNA in a Next, we focused on the secondary structure of the two DCL1 sequence-independent manner preferring flexible single-stranded cleavage sites flanking the miRNA/miRNA* duplex, which and non-Watson–Crick base-paired regions18,19. Probing the two correspond to positions −1, 1, 2, 3 and 21, 22, 23, 24 for the precursors with Tb (III) showed that the C-G variant reduced the first and the second cleavage site (Fig. 1d, e, Supplementary flexibility of the precursor (Fig. 2j). Most conspicuously, the Fig. 1) (Source Data). Watson–Crick base pairs were the most analysis revealed that three phosphodiester bonds in the vicinity frequent interactions, especially at positions 1, 3, and 23 (Fig. 1d, of the mutation were selectively stabilized in the MIR172A C-G e). In contrast to these positions, position −1 was unpaired in mutant precursor (Fig. 2j), suggesting a further stabilization of ~40% of the precursors (Fig. 1d, e). The preferred mismatch at the base pairs near the mutation. Interestingly, MIR172A this position was an A-C/C-A with an overall frequency of 25% has been shown to be the most active member of the miR172 (Arabidopsis) and 24% (eudicots) (Fig. 1d, e). Interestingly, C-C family in Arabidopsis16,17. Our results show that a single change mismatches, which were infrequent in the overall precursor can increase its processing efficiency and activity even further structures (Fig. 1b, c), were absent from DCL1 cleavage sites in vivo. (Fig. 1d, e). These biases were not explained by the frequencies of the individual nucleotides, as they were similar to each other (Supplementary Table 1). The nucleotide identity at the DCL1 cleavage site significantly As a comparison, we analyzed the composition of the cleavage affects MIR172 processing. Next, we systematically analyzed the sites of human miRNA precursors. All the positions correspond- importance of the base pair composition at position 1 of the ing to the first cleavage site were rather similar in human MIR172A pri-miRNA by generating mutant variants with the 16 precursors with ~70% of Watson–Crick base pairs (Supplemen- different possible combinations of nucleotides (Fig. 3). These tary Fig. 2). No obvious bias in the preference of the nucleotide precursors were cloned under the control of the 35 S viral pro- pairs was observed, and C-C mismatches were present at similar moter and transformed into plants. We first compared the wild- frequencies as other mismatch variants (Supplementary Fig. 2), type precursors with mutant variants harboring closed pairs, suggesting that the features we found were plant specific. either with canonical Watson–Crick base pairs or G-U wobble (Fig. 3a). All these variants accumulated miRNA levels similar to the wild-type precursor, with the exception of the MIR172A C-G A single-nucleotide polymorphism increases the processing variant that accumulated the highest levels of the small RNA efficiency of MIR172A. Of the different nucleotide pairs forming (Fig. 3a and Supplementary Fig. 5a). We also determined mature the first DCL1 cleavage site, the G-C/C-G pair at position 1 had miR172 and pri-miR172 by RT-qPCR in these mutants and the highest frequency, being present at 59% and 55% of the confirmed that the MIR172 C-G variant had the highest ratio precursors of Arabidopsis and eudicots, respectively (Fig. 1d, e). miR172/pri-miR172 (Supplementary Fig. 5b–d). Furthermore, this position is paired in 90% or more of the plant Differently to the results obtained with MIR172A mutants that miRNA precursors. Exceptions to this rule are the Arabidopsis harbored closed pairs, we found that precursors with mismatches MIR827 and MIR172A precursors (Supplementary Fig. 3). We accumulate largely variable amounts of mature miRNA levels, analyzed the importance of the pairing at position 1 in MIR172A, which can decrease >10-fold depending on the identity of which causes obvious phenotypes after overexpression, namely nucleotides forming the mismatch (Fig. 3b and Supplementary early flowering time and flower-patterning defects, that quanti- Fig. 6). Interestingly, the A-C and C-C mismatch variants, which tatively correlate with the mature miRNA levels10,11,16,17. generate miRNAs of identical sequences, accumulated largely MIR172A harbors a C-U mismatch at position 1 (Fig. 2a, Source different levels of miRNAs (Fig. 3b). Although the A-C mutant Data), so we introduced a single change by site-directed muta- accumulated nearly half the levels of mature miRNA with respect genesis in the MIR172A precursor and generated a C-G variant. to the wild-type precursor, the levels of miR172 were barely We analyzed miR172 levels in inflorescences of primary detected in the C-C mismatch variant (Fig. 3b and Supplementary transgenic (T1) plants overexpressing the wild-type (wt) and Figure 6). Furthermore, the mutant variants A-A, A-G, and G-A the C-G variant precursor from the 35S viral promoter using accumulated ~10% mature miRNA levels with respect to the C-U small RNA blots (Fig. 2b). Overexpression of the MIR172 C-G wt MIR172A (Fig. 3b). mutant variant precursor resulted in >60% increase in the We determined the levels of miR172* (Fig. 3b, lower panel) accumulation of mature miRNA with respect to control and found a similar pattern of accumulation for both guide and inflorescences transformed with wt MIR172 precursor (Fig. 2b). passenger miR172 (Fig. 3b, upper panel), suggesting that Although both precursors caused a strong acceleration of differences in miR172 levels were not caused by preferentially flowering (Supplementary Fig. 4a), we found that the C-G variant loading miR172* into an AGO complex. We analyzed the pri- produced stronger flower-patterning defects (Fig. 2c), as expected miRNA levels by RT-qPCR of the mutants that failed to from the higher levels of mature miRNA. To further validate accumulate miR172 and found that the mutant precursors these results, we determined the levels of pri-miR172 and mature accumulate larger amounts of primary transcripts (Fig. 3c). These miR172 by RT-qPCR in seedlings (Fig. 2d, f) and inflorescences results indicate that these mutants were affected at the level of the (Fig. 2g, i). Overall, we observed that the MIR172 C-G mutant precursor processing. In good agreement with the low miR172 variant accumulated higher levels of miR172 than the wt MIR172 levels, overexpression of the MIR172 C-C variant caused only a C-U precursor (Fig. 2e, h), whereas the levels of pri-miR172 minor effect on flowering time and no defects on the flower were similar or lower (Fig. 2d, g). In consequence, the ratio patterning (Fig. 3d, e). NATURE COMMUNICATIONS | (2020)11:5320 | https://doi.org/10.1038/s41467-020-19129-6 | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-19129-6 a b c j Prom35S: Flower defects 50 100 Tb(III) MIR172A Empty wt vector - – 1 1 T1 A – + C-U ++ + C-G C-U (d) (d)OH Prom35S: MIR172A – C-G + ++ ++ Replicate wt Empty 1 1 Empty vector C-U C-G vector C-U C-G 1.0 5.1 8.6 1.0 2.6 5.7 miR172 U6 EB d 1.6 e 380 c f 780 *** pri-miR172a pri-miR172a relative expression miR172 relative expression 1C-U wt Ratio miRNA/pri-miRNA ** 1C-G 0.8 190 390 b n.d a U 0.0 0 0 C Empty C-U C-G Empty C-U C-G C-U C-G G vector Prom35S: vector Prom35S: Prom35S: MIR172A MIR172A MIR172A g 2.6 h18 b i 9.0 ** pri-miR172a relative expression miR172 relative expression Ratio miRNA/pri-miRNA 1.3 9 4.5 a a n.d 0.0 0 0.0 Empty C-U C-G Empty C-U C-G C-U C-G 2 4 6 8 10 12 vector Prom35S: vector Prom35S: Prom35S: Normalized intensity (x10–4) MIR172A MIR172A MIR172A Fig. 2 A single-nucleotide modification enhances miRNA biogenesis in vivo. a Predicted pri-miR172a foldback structure. The miRNA is shown in pink and the miRNA* in light purple. Letters on the left correspond to the nucleotides indicated with arrows in j. b Small RNA blots of primary transgenic inflorescences overexpressing MIR172A wt (1C-U) and mutated, 1C-G. Numbers above the gel indicate signal intensity normalized to U6 signal and relative to plants expressing the empty vector. The ethidium bromide (EB) staining of each gel is shown at the bottom. Blots were cropped to remove irrelevant lanes; see Supplementary Fig. 11a for uncropped images of the gels. c Phenotypes of plants overexpressing MIR172A wt (1C-U) and mutated (1C-G). Percentage of plants with no (−), moderate (+), and strong (++) flower defects. A representative flower for each phenotype is indicated in the inset. d–i Accumulation of pri-miRNA (d, g) and miRNA (e, h) as well as ratios between miRNA and pri-miRNA (f, i) in seedlings (d–f) or inflorescences (g–i) expressing the empty vector, wt MIR172A (1C-U), or mutated 1C-G as determined by RT-qPCR. Each sample indicated by different shades of pink corresponds to a pool of at least ten independent transgenic plants. Asterisks in d, f, and i indicate statistically significant differences according to Student t test (two-tailed), p < 0.01 (**) and p < 0.001 (***). Different letters in e and h indicate statistically significant differences, according to ANOVA followed by Tukey’s multiple comparison test (p < 0.05). n.d., not detected. j Denaturing polyacrylamide gel of Terbium induced autohydrolysis of pri-miR172a wt (1C- U) and mutated (1C-G). From right to left: in vitro transcribed pri-miR172a incubated with Tb (III), C-G and wt (C-U) respectively; T1 (d)T1 RNase in denaturing conditions; A(d): RNAse A in denaturing conditions; OH−: alkaline hydrolysis. Two independent reactions were loaded for each treatment. The intensity profile is depicted on the left. Arrowheads indicate bands with different intensity in wt and the 1C-G mutant variant. See Supplementary Fig. 12 for uncropped image of the gel. C-C mismatches systematically affect plant miRNA processing. the importance of the secondary structure in MIRNA proces- As the identity of the unpaired nucleotides at position 1 of sing9–11. However, the C-C mismatch at position 1 reduced even MIR172A significantly affected miRNA biogenesis, we decided to further the biogenesis of miR165a than the U-C mismatch test their importance in other precursors and positions (Fig. 4). (Fig. 4a and Supplementary Fig. 7a, b). We also generated two At position −1, MIR165A has a mismatch followed by a new mutants in MIR172A, an A-C and C-C mismatch at position Watson–Crick G-C pair at position 1 (Fig. 4a). Introducing a −1. Although both mismatches affected miR172 biogenesis, the mutation at position 1 that generated two continuous mismatches C-C mismatch had a stronger effect than the A-C variant (Fig. 4c affected the biogenesis of miR165a (Fig. 4a, b and Supplementary and Supplementary Fig. 7e, f). These results confirm the dele- Fig. 7a, b), in good agreement with previous work highlighting terious impact of C-C mismatches in the first DCL1 cleavage site. 4 NATURE COMMUNICATIONS | (2020)11:5320 | https://doi.org/10.1038/s41467-020-19129-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-19129-6 ARTICLE a Prom35S:MIR172A wt Empty 1 1 1 1 1 1 1 Position vector C-U G-U U-G A-U U-A C-G G-C nt id n.d. 1.0 0.9 1.2 1.2 1.4 2.3 1.3 Position 3 2 miR172 1 -1 U6 EB b Prom35S:MIR172A Prom35S:MIR172A wt wt Empty 1 1 1 1 1 1 1 Empty 1 1 1 1 Position vector C-U U-U U-C C-C A-A A-C C-A vector C-U G-G A-G G-A nt id n.d. 1.0 0.3 0.3 <0.1 <0.1 0.4 0.1 n.d. 1.0 1.2 0.2 0.1 miR172 n.d. 1.0 0.4 0.6 0.2 0.1 1.0 0.2 n.d. 1.0 0.8 0.1 0.3 miR172* U6 EB c d a e 100 22 b pri-miR172a relative expression 15 Number of leaves Trangenics (%) ab c ab 10 50 1 ab 5 b a n.d 0 – + ++ – + ++ – + ++ 0 Empty wt 1 1 1 1 Empty 1 1 Empty 1 1 vector C-U C-C G-A A-G A-A vector C-U C-C vector C-U C-C Prom35S:MIR172A Prom35S: Prom35S: MIR172A MIR172A Next, we analyzed mutations in the second DCL1 cleavage A-C, U-C, G-U, C-C, and G-G (Fig. 4g, Supplementary Fig. 7g). site. The MIR397A precursor has an A-C mismatch at position We observed a broad distribution of miR164 accumulation 21 (Fig. 4e). We introduced C-C and A-G mismatches and depending on the identity of nucleotides (Fig. 4g, Supplementary observed that the C-C was detrimental for miR397 Fig. 7g), similar to what we observed previously for miR172 biogenesis while the A-G variant was not (Fig. 4e and (Fig. 3). We observed that C-C and G-G mutants accumulated Supplementary Fig. 7c, d). MIR164C has a G-C pair at position much less miR164 than A-C, U-C, and G-U variants (Fig. 4g and 23, and we generated five variants by site-directed mutagenesis, Supplementary Figure 7g). NATURE COMMUNICATIONS | (2020)11:5320 | https://doi.org/10.1038/s41467-020-19129-6 | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-19129-6 Fig. 3 Nucleotide pair variants control miR172 biogenesis. a, b MIR172A was modified for all the nucleotide pair variants in position 1 of the precursor. Seedlings expressing all the possible pair variants in position 1 of the precursor were analyzed by small RNA blot for miR172. Numbers above gels indicate signal intensity normalized to U6 signal and relative to the wt precursor signal (1C-U). On a, a schematic representation of MIR172A is depicted, showing the position of the base pair modification. In a, b above each gel, a portion of the precursor is included, with the specific nucleotides modified indicated with a yellow circle. Bottom panel in b, a small RNA blot for miR172a*. The ethidium bromide (EB) staining of each gel is shown at the bottom. Blot b was cropped to remove irrelevant lanes; see Supplementary Fig. 11b–d for uncropped images of the gels. See Supplementary Figs. 5 and 6 for biological replicates. c Pri-miRNA quantification by RT-qPCR of seedlings expressing MIR172A unpaired variants: wt, 1C-C, 1 G-A, 1 A-G and 1 A-A. Different letters indicate statistically significant differences (ANOVA followed by Tukey’s multiple comparison test (p < 0.05)). n.d.: not detected. d, e Phenotypes of plants overexpressing MIR172A wt (1C-U) and mutated (1C-C) variant. Note the small effect of the mutated MIR172A C-C on flowering time (d) and flower- patterning defects (e) (see Fig. 2c for phenotypes references). Different letters in d indicate statistically significant differences (Kruskal–Wallis multiple comparison test, p < 0.05). As expected, we observed a correlation between the miRNA Data 1). To check the stability of the MIRNAs through the levels generated by the different mutant precursors and the simulations, the total axial bend of each model was calculated, strength of the developmental defects (miR165, miR164) or the which showed a similar stable behavior in all cases (Supple- target transcript levels (miR397) (Supplementary Fig. 8). In mentary Data 1). These simulations were used to analyze in detail addition, all the MIRNA mutants that failed to accumulate all the structural features of all the base pairs belonging to each substantial levels of the mature miRNA expressed primary nucleic acid chain (Supplementary Data 1 and 3). The structural transcripts at similar or higher levels than their cognate wild- parameter that showed the most distinct behaviors among the type MIRNAs (Fig. 4b, d, f, h), demonstrating that the low levels different mutant variants was the shear parameter, which repre- of the mature miRNAs were not caused by a poor expression of sents the relative displacement between the bases on a plane their precursors. perpendicular to the RNA axis (Fig. 5a–d, Supplementary Overall, the results show that the identity of unpaired Data 1). nucleotides in the miRNA precursors has a strong effect on the Although canonical Watson–Crick base pairs as well as G-U accumulation of the mature miRNA in vivo (Figs. 3 and 4). wobbles showed a stable value of the shear parameter along the Although the final impact caused by the unpaired nucleotides is MD simulation (Supplementary Data 1), mismatch variants, likely affected by the precursor context, we found that C-C displayed higher fluctuation of values. Some variants, such as C-U mismatches consistently impaired miRNA biogenesis. Contrarily, and A-C showed some fluctuation between stable values. Others, the biogenesis of precursors harboring A-C, U-U, and U-C such as C-C had a higher fluctuation and no stable values at all mismatches is similar to those with WC pairs. Interestingly, the (Fig. 5d). To better illustrate the different behavior of the latter mismatches were the most frequent in miRNA precursors nucleotide pair variants, we performed a clustering of the (Fig. 1). different conformations visited through the MD simulation Next, we studied the effect of the identities of mismatches in (Fig. 5b, Supplementary Data 3) and an analysis of the formation the miRNA/miRNA* duplex region. Introducing a mismatch at and permanence of hydrogen bonds between the bases (Fig. 5c, position 6 of miR172/miR172* to generate an A-C mismatch had Supplementary Data 3). a small effect on miR172 accumulation (Fig. 4i–j, Supplementary Canonical Watson–Crick base pairs and G-U wobbles had a Fig. 7h–i and Supplementary Fig. 8d, e). We observed, however, a single stable conformation, in agreement with the shear values. different effect after introducing a C-C mismatch variant at the Variants such as G-G and A-C showed different positions for the same position. On the one hand it diminished miR172 bases along time (Supplementary Data 3), reflecting periodical accumulation (Fig. 4i and Supplementary Fig. 7h, i). On the changes between different conformations that are stabilized other hand, it caused the generation of two small RNA variants through hydrogen bonds (Supplementary Data 3). In contrast, from the MIR172A precursor (Fig. 4i). So, although C-C analysis of the C-C mismatch showed constant movement mismatches impair the processing quantitatively in all the cases between the different conformations, consistent with the lack of analyzed so far, the mutation in the miRNA/miRNA* region any stable conformation for this variant (Fig. 5b–d, Supplemen- affected the biogenesis qualitatively. It has been recently shown tary Data 1 and 3). Results for the C-C variant suggest that the flexibility on miR168/miR168* interaction can be linked significantly higher flexibility than in the wild-type precursor to miR168 variants of different lengths and sequences20. We harboring a C-U mismatch. We also analyzed the variants at inspected the miR168/miR168* duplex and found a C-C position 23 of MIR164C (Fig. 4g and Supplementary Fig. 10a) by mismatch that is conserved in miR168/miR168* duplexes in 37 MD. In good agreement, the shear parameter showed again an out of 38 precursors (Supplementary Fig. 9), so that conservation unstable behavior of the C-C pair (Supplementary Fig. 10b and of this particular mismatch is also consistent with particular Supplementary Data 2 and 3). properties of the precursor processing. In order to confirm these results, we used Terbium (III) to probe the structure of the wild-type MIR172A Molecular dynamics analysis revealed differences in the mis- precursor harboring C-U and the mutant with a C-C mismatch. match variants of MIR172A. The surprisingly different effects on Despite both precursors having a mismatch at position 1, we miRNA biogenesis caused by the identity of nucleotides promp- found that the MIR172A C-C precursor variant was more ted us to study the system in more detail. We resorted to mole- sensitive to Terbium (Tb III) than the wild-type precursor at cular dynamics (MD) simulations in order to obtain atomistic the location of the C-C mismatch as well as at neighboring details about the influence of the base pair nature on the pre- nucleotides (Fig. 5e), supporting our MD calculations, and cursor behavior. To achieve this goal, the 16 variants for position demonstrating an increase of structural flexibility caused by the 1 of a segment of the MIR172A precursor (Figs. 3 and 5a) were C-C mismatch. built in silico and subjected to 100 ns of MD simulations. A In addition, we mapped the processing intermediates for representative snapshot of the system was built for the wt variant, MIR172A with a C-C mismatch and observed a random composed of the dsRNA, solvent box and NaCl (Supplementary distribution of the cuts (Fig. 5f), in contrast to the wild-type 6 NATURE COMMUNICATIONS | (2020)11:5320 | https://doi.org/10.1038/s41467-020-19129-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-19129-6 ARTICLE a c e Prom35S:MIR165A Prom35S:MIR172A Prom35S:MIR397A Position 24 23 22 wt wt 21 wt Empty 1 1 1 Position Empty –1 –1 –1 Position Empty 21 21 21 Position vector G-C U-C C-C nt id vector G-C A-C C-C nt id vector A-C C-C A-G nt id 1 16 2 1 n.d. 1.0 0.2 <0.1 n.d. 1.0 0.2 0.8 3 3 2 miR165 2 miR172 miR397 1 1 –1 –1 Position Position U6 U6 U6 EB EB EB b 150 b d 6.0 f 2.0 pri-miR397a relative expression pri-miR165a relative expression pri-miR172a relative expression b a b b ab b 75 3.0 1.0 b a a n.d n.d 0 0.0 0.0 Empty G-C U-C C-C Empty G-C A-C C-C Empty A-C C-C A-G vector vector vector Prom35S:MIR165A Prom35S:MIR172A Prom35S:MIR397A g i Prom35S:MIR164C Prom35S:MIR172A Position 24 23 22 21 Empty Empty wt vector wt Position vector 23 23 23 23 23 23 6 6 6 Position G-C A-C U-C G-U C-C G-G nt id G-C C-C A-C nt id n.d. 1.0 0.3 0.6 0.3 <0.1 <0.1 7 1 30 1 10 6 miR164 5 Position miR172 U6 U6 EB h 9.0 j 5.0 pri-miR164c relative expression pri-miR172a relative expression 4.5 2.5 n.d n.d 0.0 0.0 Empty G-C A-C U-C G-U C-C G-G Empty G-C C-C A-C vector vector Prom35S:MIR164C Prom35S:MIR397A Fig. 4 miRNA biogenesis is strongly affected by the nucleotide identity at the precursor unpaired positions. Several miRNA precursors where modified and overexpressed in plants. Schemes of miRNA precursors are included to the left in each panel, with the miRNA in pink and miRNA* in light purple. Seedlings expressing the MIRNA variants were analyzed by small RNA blots (a, c, e, g, i) and RT-qPCR (b, d, f, h, j). Numbers above gels indicate signal intensity normalized to U6 signal and relative to empty vector (a, i) or wt precursor on each gel (c, e, g). The ethidium bromide (EB) staining of each gel is shown below the blots. See Supplementary Fig. 7 for biological replicates and Supplementary Fig. 11e–i for uncropped gel images. The pri-miRNA levels were determined in three biological replicates by RT-qPCR and indicated below the blots. Different letters indicate significant differences (ANOVA followed by Tukey’s multiple comparison test (p < 0.05)). n.d., not detected. a, b MIR165A wt (1 G-C) and 1 U-C and 1 C-C mutant variants. c, d MIR172A wt (−1 G-C), and -1 A-C and −1 C-C mutant variants. e, f MIR397A wt (21 A-C) and 21 C-C and 21 A-G mutant variants. g, h MIR164C wt (23 G-C) and 23 A-C, 23 U-C, 23 G-U, 23 C-C, and 23 G-G mutant variants. i, j MIR172A wt (6 G-C) and 6 C-C and 6 A-C mutant variants. NATURE COMMUNICATIONS | (2020)11:5320 | https://doi.org/10.1038/s41467-020-19129-6 | www.nature.com/naturecommunications 7 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-19129-6 a b e Tb(III) wt MIR172A C-U 1 1 T1 A - C:3 (d) (d) OH C-U C-C –4 Normalized intensity (×10 ) Cluster C:2 2 4 6 8 10 12 C:1 C:0 Shear 0 20 40 60 80 100 x Time (ns) U C G MIR172A C-C U y C:3 z Cluster C:2 C:1 Sequence used for MD simulation C:0 0 20 40 60 80 100 Time (ns) c MIR172A C-U MIR172A C-C Cluster Cluster 1.96% a 54.4% b a C:0 4.68% C:0 pri-miR172a 1.16% d 1C-U wt e 4.32% 1C-C C:1 C:1 C:2 U C:2 C G U C:3 d MIR172A C-U 8 4 Shear (Å) 0 –4 20 40 60 80 100 20 40 60 80 100 0 20 40 60 80 100 Time (ns) Time (ns) Time (ns) MIR172A C-C f 8 MIR172A 5´ 4 4/12 Shear (Å) 0 wt 1C-U –4 3´ 8/12 –8 0 20 40 60 80 100 0 20 40 60 80 100 0 20 40 60 80 100 5´ 1/10 1/10 Time (ns) Time (ns) Time (ns) 1/10 1C-C 1/10 2/10 3´ 2/10 1/10 1/10 Fig. 5 A C-C variant increases the flexibility of the precursor structure. a Schematic representation of MIR172A precursor. Letters on the left correspond to the nucleotides (nt) indicated with arrows in e. A dashed outline indicates the region used for molecular dynamics (MD) simulations, the sequence is showed below, “x” indicates the nucleotides modified for the simulation of each variant. The “shear” parameter (an in-plane shift of the bases within the base pair) is depicted on the right (movement of bases along the x axis). See Methods and Supplementary Data 1 for additional information. b–d Results of MD for MIR172A 1C-U (wt) and 1C-C variants in position 1. Clustering (b), nucleotide pair conformation (c), and shear parameter (d) are depicted. The percentage of simulation time that each hydrogen bond is formed is indicated next to the bonds in the conformational depictions (c), numbers indicating the conformations correspond to clusters in b. For shear parameter three replicates are shown. See Supplementary Data 1 for all other nucleotide variants. e Structural analysis of pri-miR172a with wt (C-U) and C-C variant at position 1. Denaturing polyacrylamide gel, from left to right: T1(d): T1 RNase in denaturing conditions; A(d): RNAse A in denaturing conditions; OH−: alkaline hydrolysis; pri-miR172a incubated with Tb(III), wt and C-C variants respectively; two independent reactions were loaded for each treatment. On the right, the intensity profile of Terbium reactions is depicted. Black arrowheads indicate bands with different intensity between wt and C-C variant. See Supplementary Fig. 12 for uncropped image of the gel. f MIR172A wt and 1C-C variant where analyzed through modified 5´RACE. Arrowheads and numbers indicate processing intermediates identified by sequencing and their frequency; green arrowheads indicate the expected DCL1 cleavage product. 8 NATURE COMMUNICATIONS | (2020)11:5320 | https://doi.org/10.1038/s41467-020-19129-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-19129-6 ARTICLE precursor that harbored ~70% of the cuts flanking the miRNA/ We found that positions 1, 3, and 23 of the precursors are miRNA*. Altogether, the results show that the identity of the mostly paired, and showed that closing the secondary structure of unpaired nucleotides confer different structural constraints that MIR172A that has an unusual mismatch at position 1 increased can be translated into an aberrant or altered miRNA precursor miR172 accumulation. In a similar way, introducing a mismatch processing. at position 23 of MIR164C systematically reduced its biogenesis. In good agreement, recent results have shown that the processing of the precursors is also favored by a G-C-rich signature in the Discussion miRNA/miRNA* duplex31, highlighting the importance of a MiRNAs are processed from precursors harboring imperfect stable secondary structure at certain regions of the plant miRNA foldback structures. Many reports have demonstrated the precursors. importance of the precursor secondary structure during miRNA Structure–function analyses of RNA foldbacks usually focus on biogenesis5–12. These structural determinants are mostly repre- the secondary structure, disregarding the importance of the sented by dsRNA segments of 15–17 bp located below or above nucleotide identity at unpaired positions. Artificial miRNA the miRNA/miRNA*5–12. In addition, the identity of the 5′ (amiRNA) technology is used to downregulate gene expression, nucleotide as well as mismatches in the miRNA/miRNA* region and is based on the replacement of the small RNA sequence of a have been shown to be important for selection of the specific natural miRNA precursor maintaining its secondary structure32– AGO protein21–26. 34. amiRNAs are especially suited to downregulate gene expres- Here, we describe a new set of sequence features that are sion in specific tissues by driving their expression from different important to define a plant miRNA primary transcript. We promoters. Our findings show the importance of considering the identified certain preferred nucleotide pairs at the precursors’ identity of the unpaired nucleotides during the design of these cleavage site, which might contribute to the recognition and/or experiments. In the case of artificial miRNAs, we suggest activity of the processing complex. However, an unexpected new designing precursors harboring only A-C, U-C, or U-U mis- finding is that the identity of the nucleotides located at mis- matches. Our results could be useful for the design of targeted matched positions of the imperfect dsRNA that serves as a modifications of MIRNA genes in order to modify the expression miRNA precursor can have a deep impact on the processing levels without altering the small RNA sequence. Furthermore, it efficiency. This is in turn reflected in the evolution of MIRNAs, as might be interesting to test whether the identity of the mismatch certain mismatches such as C-C and G-G are largely excluded pairs, especially C-C mismatches, also affects other biological from the precursor structures, especially from DCL1 cleavage processes. sites. Although we did not find a similar bias in the nucleotide pairs at the cleavage sites of human precursors, a general analysis Methods of the frequency of nucleotide pairs in other RNA structures Sequence analysis of MIRNAs. Analysis of nucleotide pairs in the precursors’ shows that the C-C mismatch is the least frequent pair27–29, sequences was performed for MIRNAs conserved in eudicots, using publicly suggesting that it might have an adverse effect in other available sequence data12. Sixty-one sequences of conserved precursors from Arabidopsis thaliana, and 1681 sequences of precursors from eudicot species were pathways too. analyzed (Supplementary Fig. 1 and Source Data). Secondary structures of miRNA We found that, at least during plant miRNA biogenesis, the precursors were predicted with the Mfold web server (http://unafold.rna.albany. C-C pair is detrimental in different locations within a miRNA edu/?q=mfold/RNA-Folding-Form) using default parameters35. For human precursor. A MD analysis showed that bases at this pair are in miRNA analysis, information of 477 precursors available in miRBase v22 were included36. Data used for these analyses is available in Source Data file. continuous movement (shear), while Terbium (III) probing revealed that a dsRNA containing a C-C mismatch is more flexible than one harboring a C-U mismatch. Taken together, Plant material. A. thaliana ecotype Columbia was used for all experiments. Plants were grown at 22 °C under continuous light at 100 µmol photons m−2 s−1. these results show that mismatches formed by different nucleo- Transgenic plants were selected on Murashige and Skoog media with 50 µg/ml tides can affect the structural properties of the foldbacks, which in kanamycin. For the experiments included in Figs. 2c, 3d–e, and Supplementary turn is reflected in a differential miRNA processing efficiency. Figs 4a, 8d, e, transgenic plants were transferred to soil after seven days and Other mismatches such as G-G, A-G/G-A are also found at low continued growing in long days. Upon flowering, primary inflorescences were frequency in plant miRNA precursors and, at least in certain collected and phenotypic analysis was performed (see Source Data file for addi- tional details). Seedlings were collected after 10 days for RNA extraction (experi- contexts, also impaired miRNA biogenesis. Mutations in the first ments included in Figs. 3–5 and related Supplementary Figures). cleavage site of the MIR172A precursor that impaired miRNA biogenesis resulted in an increased accumulation of the pri- Transgenic lines. miRNA precursors corresponding to MIR165A (AT1G01183), miR172 levels, indicating that the first cleavage reaction was MIR164C (AT5G27807), MIR172A (AT2G28056), and MIR397A (AT4G05105) compromised, however we cannot rule out the C-C mismatches were expressed under the control of the CaMV 35 S promoter of the binary vector interfere with other steps of the miRNA pathway as well, such as CHF337. Sequences of all the constructs used and prepared in this article are AGO loading. indicated in Supplementary Table 3. Although C-C mismatches are usually excluded from the precursors’ foldbacks, we found a few cases where C-C mis- RNA analysis. Ten-day old seedlings or inflorescences were ground in liquid matches are evolutionarily conserved in the precursor structure, nitrogen, and RNA extraction was performed with TRI Reagent® (MRC) following the product instructions. Each sample analyzed contained material from 10 inde- such as the one found in the miRNA/miRNA* duplexes of pendent primary transgenic plants. Small RNA blots were performed as described MIR168. The later precursor has been shown to produce isomiRs previously38. Probes used for detection of each miRNA are listed in Supplementary with different features, which in turn is relevant for the selection Table 3. Note that we used 19 or 20 nt long probes that perfectly match to all the of the cognate AGO protein20,30. A flexible miR168/miR168* mutant miRNAs generated in our experiments. In brief, 3–15 µg of total RNA were resolved on 17% polyacrylamide denaturing gels (7 M Urea), transferred to a duplex structure has been shown to be responsible for the miR168 Nytran SPC membrane (GE Healthcare) and UV-crosslinked. Probes were 5′-end- isomiRs20. In this line, the flexibility introduced by a C-C mis- labeled with [γ-32P] ATP using T4 polynucleotide kinase (Thermo Scientific). match to dsRNA segments might fulfill specific roles in certain Detection of miR172a* was carried out first, the membraned stripped and checked, specific contexts. Of all the MIR172A variants tested, the only one and then used for miR172. Images were obtained with a Typhoon FLA 7000 (GE Healthcare) and the signal intensity of each band was determined using the Gel- that generated isomiRs was a precursor harboring a C-C mis- Quant.NET software (biochemlabsolutions.com). The intensity values of the bands match in the miRNA/miRNA* duplex, thus supporting this detected with the miRNA probes were normalized to the signal obtained with the possibility. U6 snRNA probe. These values were then divided by the normalized signal NATURE COMMUNICATIONS | (2020)11:5320 | https://doi.org/10.1038/s41467-020-19129-6 | www.nature.com/naturecommunications 9 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-19129-6 intensity of the sample corresponding to the transgenics with the empty vector or Reporting summary. Further information on research design is available in the Nature with the wt precursor. All the quantifications are detailed in Source Data file. Research Reporting Summary linked to this article. Primary transcripts and small RNAs were measured by RT-qPCR in at least three independent samples (two technical replicates) as described before12. In brief, 500 ng of total RNA were treated with RQ1 RNAse-free DNase (Promega). First Data availability strand synthesis was carried out with M-MLV reverse transcriptase (Invitrogen) All relevant data are included in the paper and/or its Supplementay Information using dTV oligonucleotides for polyadenylated RNAs and stem-loop files. Data corresponding to Figures 1–5, Supplementary Figure 2 and Supplementary oligonucleotides for mature miRNAs or miR172a*39,40. PCR reactions were Figures 4-9 are detailed in Source Data file provided with this manuscript. performed in an AriaMx Real-time PCR System (Agilent) using SYBR Green I (Roche) to monitor double-stranded DNA synthesis. Pri-miRNA levels were Received: 23 January 2020; Accepted: 24 September 2020; normalized to RPS26 (AT3G56340) expression. Oligonucleotides used are listed in Supplementary Table 2. The rapid amplification of 5′ cDNA ends (5′ RACE) method to detect cuts in the precursors was carried out as described previously14. In brief, total RNA was treated with DNAse, purified and the 5′-end of processing intermediates were ligated to a RNA adapter (GeneRacer™ RNA, Thermofisher). Then, polyadenylated RNA was purified with PolyATtract® mRNA Isolation Systems (Promega). First References strand synthesis was carried out with an oligonucleotide specific for transcripts 1. Yu, Y., Jia, T. & Chen, X. The ‘how’ and ‘where’ of plant microRNAs. N. derived from CHF3 vector. Hemi-nested PCRs were performed with two vector- Phytol. 216, 1002–1017 (2017). specific reverse oligonucleotides. PCR products were purified, cloned, and 2. Achkar, N. P., Cambiagno, D. A. & Manavella, P. A. miRNA biogenesis: a sequenced. Oligonucleotides used are listed in Supplementary Table 2. dynamic pathway. Trends Plant Sci. 21, 1034–1044 (2016). 3. Ma, Z. & Zhang, X. Actions of plant Argonautes: predictable or unpredictable? Curr. Opin. Plant Biol. 45, 59–67 (2018). 4. Fang, X. & Qi, Y. RNAi in plants: an Argonaute-centered view. Plant Cell 28, RNA structure. Pri-miR172a wt (C-U), 1C-G, and 1C-C versions were transcribed 272–285 (2016). in vitro using T7 RNA polymerase as described before6. Alkaline hydrolysis of 5. Bologna, N. G. et al. Multiple RNA recognition patterns during microRNA transcripts was used as a ruler, as well as MW standards. Terbium (Tb III) induced biogenesis in plants. Genome Res. 23, 1675–1689 (2013). autohydrolysis reactions were carried out by incubating 1.5 µl RNA (1 µM), 1 µl folding buffer (20 mM Tris pH 7.5, 25 mM KCl and 10 mM MgCl2), 0.3 µl of 6. Moro, B. et al. Efficiency and precision of microRNA biogenesis modes in Escherichia coli tRNA (0.1 mg/ml) and 6.2 µl water for 15 min at room tempera- plants. Nucleic Acids Res. 46, 10709–10723 (2018). ture, after which 1 µl of 100 mM Tb (III) was added. The reaction was carried out 7. Zhu, H. et al. Bidirectional processing of pri-miRNAs with branched terminal for 15 min, and then stopped by precipitating the RNA overnight. For control, 1 µl loops by Arabidopsis Dicer-like1. Nat. Struct. Mol. Biol. 20, 1106–1115 (2013). of water was added instead of Tb (III). After RNA precipitation and resuspension, 8. Cuperus, J. T. et al. Identification of MIR390a precursor processing-defective the digested products were analyzed in 8% or 12% (w/v) denaturing polyacrylamide mutants in Arabidopsis by direct genome sequencing. Proc. Natl Acad. Sci. gels, respectively. The results were visualized by Phosphor Imaging (Typhoon FLA USA 107, 466–471 (2010). 7000, GE). Image processing, band profile acquisition and quantification were 9. Song, L., Axtell, M. J. & Fedoroff, N. V. RNA secondary structural performed with the FIJI software41. determinants of miRNA precursor processing in Arabidopsis. Curr. Biol. 20, 37–41 (2010). 10. Werner, S., Wollmann, H., Schneeberger, K. & Weigel, D. Structure determinants for accurate processing of miR172a in Arabidopsis thaliana. Statistical methods. Tests applied are indicated in the Figure captions. Infor- Curr. Biol. 20, 42–48 (2010). mation about size of samples (n), and the respective tables of each statistical test are 11. Mateos, J. L., Bologna, N. G., Chorostecki, U. & Palatnik, J. F. Identification of included in Source Data file. RT-qPCR and plant phenotype data were analyzed microRNA processing determinants by random mutagenesis of Arabidopsis with InfoStat software42. MIR172a precursor. Curr. Biol. 20, 49–54 (2010). 12. Chorostecki, U. et al. Evolutionary footprints reveal insights into plant MicroRNA biogenesis. Plant Cell 29, 1248–1261 (2017). MD simulations. The initial structural models for the dsRNAs were built from a 13. Addo-Quaye, C. et al. Sliced microRNA targets and precise loop-first portion of the sequences of MIR172A from position −10 to position 10 and processing of MIR319 hairpins revealed by analysis of the Physcomitrella MIR164C from position 13 to position 32 using NAB43 (Fig. 5a and Supple- patens degradome. RNA 15, 2112–2121 (2009). mentary Fig. 10a). Each model was immersed in a truncated octahedral periodic 14. Bologna, N. G., Mateos, J. L., Bresso, E. G. & Palatnik, J. F. A loop-to-base box with a minimum solute-wall distance of 8 Å filled with around 10,000 water processing mechanism underlies the biogenesis of plant microRNAs miR319 molecules and 0.1 M NaCl. Each system was then subjected to 50 ns of MD and miR159. EMBO J. 28, 3646–3656 (2009). simulation (detailed below). Subsequently, the nucleotide pair at position 1 of 15. Kurihara, Y. & Watanabe, Y. Arabidopsis micro-RNA biogenesis through MIR172A and position 23 of MIR164C were mutated in silico to generate the 16 Dicer-like 1 protein functions. Proc. Natl Acad. Sci. USA 101, 12753–12758 different possible combinations of the 4 nucleotides. All models were further (2004). subjected to 100 ns of MD simulations, during which their structures were ana- 16. Aukerman, M. J. & Sakai, H. Regulation of flowering time and floral organ lyzed. MD simulations were performed with the Amber 14 package44. The force identity by a MicroRNA and its APETALA2-like target genes. Plant Cell 15, field used was ff99bsc0 with χOL3 modifications45,40 for dsRNAs and TIP3P for 2730–2741 (2003). water46. Particle-mesh Ewald47 was implemented for long range interactions and a 17. Chen, X. A microRNA as a translational repressor of APETALA2 in cutoff distance of 12 Å was used. Temperature and pressure were regulated with Arabidopsis flower development. Science 303, 2022–2025 (2004). Berendsen thermostat and barostat48 with a coupling constant of 2 ps. In addition, 18. Choudhary, P. K., Gallo, S. & Sigel, R. K. Monitoring global structural changes we performed two extra versions of each MD simulation of MIR172A using and specific metal-ion-binding sites in RNA by in-line probing and Tb(III) Langevin thermostat with a collision frequency of 5 ps−1 to improve statistical cleavage. Methods Mol. Biol. 1086, 143–158 (2014). significance and ensure results are not dependent on thermostat choice. A com- 19. Harris, D. A. & Walter, N. G. Probing RNA structure and metal-binding sites parison of the calculated shear parameter (see below) is shown in Fig. 5d and using terbium(III) footprinting. Curr. Protoc. Nucleic Acid Chem. Chapter 6, 8 Supplementary Data 1, in which the reproducibility of the results can be seen. The (2003). simulation protocols, validated in previous work49, started with an energy mini- 20. Iki, T. et al. Structural flexibility enables alternative maturation, mization of the initial structure, followed by a heating phase from 0 K to 300 K ARGONAUTE sorting and activities of mir168, a global gene silencing through 20 ps at constant volume and an equilibration phase at 300 K and 1 bar of regulator in plants. Mol. Plant 6, 1008–1023 (2018). constant pressure through 20 ps to reach proper density. These equilibrated 21. Endo, Y., Iwakawa, H. O. & Tomari, Y. Arabidopsis ARGONAUTE7 selects structures were the starting point for 100 ns of MD simulations at 300 K in the miR390 through multiple checkpoints during RISC assembly. EMBO Rep. 14, NVT ensemble. To avoid artifacts of unpairing at the base pairs of the ends of the 652–658 (2013). models, a simple distance restraint was applied to base pairs located at the ends of 22. Zhang, X. et al. ARGONAUTE PIWI domain and microRNA duplex structure the strands so that their hydrogen bonds are maintained. Structural parameters for regulate small RNA sorting in Arabidopsis. Nat. Commun. 5, 5468 (2014). the dsRNAs were calculated using Curves+50, which include axial bend, shear, buckle, stretch, propeller, stagger, opening, shift, tilt, slide, rill, rise, and twist 23. Zhu, H. et al. Arabidopsis Argonaute10 specifically sequesters miR166/165 to (Supplementary Data 1 and 2). Base pair conformational clustering was performed regulate shoot apical meristem development. Cell 145, 242–256 (2011). using the hierarchical agglomerative approach from the cpptraj module of 24. Mi, S. et al. Sorting of small RNAs into Arabidopsis argonaute complexes is Amber14 (Supplementary Data 3). directed by the 5′ terminal nucleotide. Cell 133, 116–127 (2008). 10 NATURE COMMUNICATIONS | (2020)11:5320 | https://doi.org/10.1038/s41467-020-19129-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-19129-6 ARTICLE 25. Montgomery, T. A. et al. Specificity of ARGONAUTE7-miR390 interaction 49. Drusin, S. I., Suarez, I. P., Gauto, D. F., Rasia, R. M. & Moreno, D. M. dsRNA- and dual functionality in TAS3 trans-acting siRNA formation. Cell 133, protein interactions studied by molecular dynamics techniques. Unravelling 128–141 (2008). dsRNA recognition by DCL1. Arch. Biochem. Biophys. 596, 118–125 (2016). 26. Takeda, A., Iwasaki, S., Watanabe, T., Utsumi, M. & Watanabe, Y. The 50. Lavery, R., Moakher, M., Maddocks, J. H., Petkeviciute, D. & Zakrzewska, K. mechanism selecting the guide strand from small RNA duplexes is different Conformational analysis of nucleic acids revisited: curves+. Nucleic Acids Res. among argonaute proteins. Plant Cell Physiol. 49, 493–500 (2008). 37, 5917–5929 (2009). 27. Davis, A. R. & Znosko, B. M. Thermodynamic characterization of single mismatches found in naturally occurring RNA. Biochemistry 46, 13425–13436 (2007). Acknowledgements 28. Davis, A. R., Kirkpatrick, C. C. & Znosko, B. M. Structural characterization of We thank Daniela Liebsch and members of the Palatnik lab for comments and discus- naturally occurring RNA single mismatches. Nucleic Acids Res. 39, 1081–1094 sion. We thank Silvana Spinelli and Ramiro Rodriguez for MIR164C and MIR397A (2011). vectors. Also, we thank the technical support received from Diego Aguirre. We thank 29. Danaee, P. et al. bpRNA: large-scale automated annotation and analysis of different agencies for funding: Bunge and Born fellow (to U.C.) and CONICET fellows RNA secondary structure. Nucleic Acids Res. 46, 5381–5394 (2018). (to A.M.L.R., S.I.D., J.L.M., N.B., B.M., E.G.B.); and J.F.P., R.M.R., and D.M.M. are 30. Maunoury, N. & Vaucheret, H. AGO1 and AGO2 act redundantly in miR408- members of CONICET. Research was supported by grants from the Argentinean Min- mediated plantacyanin regulation. PLoS ONE 6, e28729 (2011). istry of Science (PICT-2015-3557) and International Centre for Genetic Engineering and 31. Narjala, A., Nair, A., Tirumalai, V., Hari Sundar, G. V. & Shivaprasad, P. V. A Biotechnology (ICGEB CRP/ARG17-01). conserved sequence signature is essential for robust plant miRNA biogenesis. Nucleic Acids Res. 48, 3103–3118 (2020). Author contributions 32. Schwab, R., Ossowski, S., Riester, M., Warthmann, N. & Weigel, D. Highly Bioinformatic analysis: A.M.L.R. and U.C. Analysis of mutant MIRNAs in Arabidopsis: B. specific gene silencing by artificial microRNAs in Arabidopsis. Plant Cell 18, M., A.S., E.G.B., N.B., J.L.M., and A.M.L.R. In vitro precursor structural studies: A.M.L.R. 1121–1133 (2006). aided by B.M. and R.M.R. MD: S.I.D. and D.M.M. Experimental design 33. Alvarez, J. P. et al. Endogenous and synthetic microRNAs stimulate A.M.L.R., S.I.D., D.M.M., and J.F.P. Wrote the manuscript: A.M.L.R., S.I.D., D.M.M., simultaneous, efficient, and localized regulation of multiple targets in diverse and J.F.P. species. Plant Cell 18, 1134–1151 (2006). 34. Carbonell, A. et al. New generation of artificial MicroRNA and synthetic trans-acting small interfering RNA vectors for efficient gene silencing in Competing interests Arabidopsis. Plant Physiol. 165, 15–29 (2014). The authors declare no competing interests. 35. Zuker, M. Mfold web server for nucleic acid folding and hybridization prediction. Nucleic Acids Res. 31, 3406–3415 (2003). 36. Kozomara, A., Birgaoanu, M. & Griffiths-Jones, S. miRBase: from microRNA Additional information Supplementary information is available for this paper at https://doi.org/10.1038/s41467- sequences to function. Nucleic Acids Res. 47, D155–D162 (2019). 020-19129-6. 37. Jarvis, P. et al. An Arabidopsis mutant defective in the plastid general protein import apparatus. Science 282, 100–103 (1998). Correspondence and requests for materials should be addressed to J.F.P. 38. Wang, L. et al. miR396-targeted AtGRF transcription factors are required for coordination of cell division and differentiation during leaf development in Peer review information Nature Communications thanks Michael Axtell, Xuhui Huang, Arabidopsis. J. Exp. Bot. 62, 761–773 (2011). and Beixin Mo for their contribution to the peer review of this work. Peer review reports 39. Chen, C. et al. Real-time quantification of microRNAs by stem-loop RT-PCR. are available. Nucleic Acids Res. 33, e179 (2005). 40. Varkonyi-Gasic, E. & Hellens, R. P. Quantitative stem-loop RT-PCR for Reprints and permission information is available at http://www.nature.com/reprints detection of microRNAs. Methods Mol. Biol. 744, 145–157 (2011). 41. Schindelin, J. et al. Fiji: an open-source platform for biological-image analysis. Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in Nat. Methods 9, 676–682 (2012). published maps and institutional affiliations. 42. Di Rienzo, J. et al. InfoStat v.2018 (Centro de Transferencia InfoStat, FCA, Universidad Nacional de Córdoba, Argentina, 2018). 43. Macke, T. J. & Case D. A. Modeling unusual nucleic acid structures. In ACS Symposium Series (ACS Publications, Washington D.C., 1998). Open Access This article is licensed under a Creative Commons 44. Case, D. A. et al. Amber 14 (University of California, San Francisco, 2014). Attribution 4.0 International License, which permits use, sharing, 45. Pérez, A. et al. Refinement of the AMBER force field for nucleic acids: adaptation, distribution and reproduction in any medium or format, as long as you give improving the description of α/γ conformers. Biophys. J. 92, 3817–3829 appropriate credit to the original author(s) and the source, provide a link to the Creative (2007). Commons license, and indicate if changes were made. The images or other third party 46. Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W. & Klein, M. material in this article are included in the article’s Creative Commons license, unless L. Comparison of simple potential functions for simulating liquid water. J. indicated otherwise in a credit line to the material. If material is not included in the Chem. Phys. 79, 926–935 (1983). article’s Creative Commons license and your intended use is not permitted by statutory 47. Brock A., Luty I. G. T. and W. F. van Gunsteren. Lattice-sum methods for regulation or exceeds the permitted use, you will need to obtain permission directly from calculating electrostatic interactions in molecular simulations. J. Chem. Phys. the copyright holder. To view a copy of this license, visit http://creativecommons.org/ 103 (1995). licenses/by/4.0/. 48. Berendsen, H. J., Jv, Postma, van Gunsteren, W. F., DiNola, A. & Haak, J. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 81, 3684–3690 (1984). © The Author(s) 2020 NATURE COMMUNICATIONS | (2020)11:5320 | https://doi.org/10.1038/s41467-020-19129-6 | www.nature.com/naturecommunications 11